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Introduction 

This report describes a method to design feedback control systems 
with both series and feedback compensators. When only the output state 
of a linear system is available, the use of both series and feedback 
condensation provides more flexibility in the design than using a 
single compensation network. The sensitivity of the closed— loop system 
to plant variations can be reduced below the values achieved by using 
only series compensation, or open loop compensation (no feedback). 
Comparisons of just series compensation with the combined series and 
feedback compensation. are made in appendices 7 and 8. 

This report begins with a description of the classical Gullleman- 
Truxal^ ^ design method and a brief description of the state variable 
feedback method of compensation More detail can be found in the 
references cited. The following section discusses the more general 
modified observer system and discusses the design procedure associated 
with modified observer systems. Several case studies and detailed 
descriptions of four numerical techniques are included as appendices. 

Three digital computer programs (FORTRAN IV) are included to 
facilitate the required manipulations. These programs are augmented by 
two supplementary programs to modify the input data and the output data. 

If the plant has N poles, then both the series compensator and 
the feedback compensator will have N-1 poles and N-1 zeros. Under these 
circumstances the plant transfer function and the desired closed loop 
transfer function are sufficient to specify all the parameters in both 
compensators except the pole positions of the feedback compensator. 

Thus, the user is free to select these pole positions to achieve the 
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lowest sensitivity. The first program leaves this choice to the user. 
The second program selects the pole positions which minimizes the 
sensitivity to low frequency variations in the plant parameters. This 
program is restricted to those plants with no zeros. The third program 
evaluates the weighted Integral sensitivity over all frequencies for a 
system with a specified compensation. The user is free to weight those 
frequencies where the disturbance is greatest. The fourth program uses 
the pole posifions of the feedback compensator as input. The output 
specifies the complete compensation and the position of the poles of 
both the series and feedback compensators . It is useful for determining 
the stability of the compensators and checking the accuracy of the 
p revious results. 


Guilleman Truxal Design Procedure 

Assume that the overall open loop gain function is given by G(s) 
as shown in Fig. 1. Then with unity feedback, the closed loop transfer 
function is given by 


Y(s) • _ G(s) 
R(s) " 1 + G(s) 


If G(s) is written as the ratio of two polynomials 


G (s) 


GN(s) 

GD(s) 


then 


GN(s) 

GD(s) + GN(s) 


Hsl 

R(s) 
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Equating the numerator and denominator polynomials 

GN(s) « C(s) 

GD(s) = R(s) - Y(s) 


( 1 ) 


The design procedure is as follows 

1. From the specifications, choose the desired closed 
loop response, Y(s)/R(s). 

2. From step 1. and eqn. 1 determine GN(a) and GD(s). 

3. Choose a suitable series conq>ensator which shifts the 
poles and zeros of the plant to the. positions specified 
by step 2. 

The above steps are simplified if one uses Gullleman’s sugges- 
tion to select a C(s)/R(s) which results in real zeros for GD(s). Then 
the zeros of GD(s) can be simply determined from a plot of R(s) - C(s) 
for only negative real values of s. One pair of complex conjugate 
poles and a group of poles on the negative real axis is usually chosen 
to meet the above requirements. 

For example suppose the plant is given by 

“ 7(fel) 

and that the specifications can be satisfied if 

100 


Ikl 

R(s) 


s^ + 20s + 100 
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From (1) 


GN(s) »= 100 

GD(s) = + 20s 


( 3 ) 


The aeries compensator C(s) and the plant transfer function must equal 
G (s) , 


C(s) Gp(s) 


100 

8 (s+20) 


From (2) 


C(s) 


100 (s+1) 
(s+20) 


(4) 


Effectively the series compensator moved the time constant of the 
motor from 1. to 1/20. 

The success of this procedure depends upon the accurate modeling 
of the motor over the range of frequencies from zero to well above 20 
radians/sec. Otherwise additional poles in the plant transfer function 
will alter the actual closed loop transfer function. In particular the 
real system may become unstable. 

State Variable Feedback Design Procedure 

This design procedure is similar in many respects to the above 
procedure. First a desired closed loop response is selected which 
satisfies the bandwidth, rise time, steady state error and other require- 
ments. The plant is divided into a set of first order transfer functions 
which may be interconnected. The output of each of these first order 
systems is fed back to the input through a fixed gain element. By 
adjusting these gain^ any desired closed loop transfer function can be 
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achieved. 

The above process is somewhat simplified if we first introduce 
a new function, Heq(s)o This function specifies the transfer function 
between the output and the input which yields the same result as state 
variable' feedback. From Fig, 2 the equivalent feedback from the output 
is 

He,(s) . '‘I + * Gj • • • 

In terms of Heq(s) the state variable feedback results in the closed 
loop response given by 


Y(s) ^ KGp(s) 

R(s) ° 1 + KGp(s) Heq(s) 


The procedure is as follows : 

1. Select the desired Y(8)/R(s) from the specifications. 

2. Equate the desired Y(s)/R(s) to the one achieved by the 
feedback from each state. 

3. Solve for the feedback coefficients. 


Details of the design procedure when zeros are present are given by 

[ 2 ] 

Schultz and Melsa . For. exanq>le if 


Gp(s) 


1 

s (s+1) 


Y(s) 

R(s) 


100 

8^ + 20s + 100 


Then if (s) » ~ , G, (s) 

X 8 Z 


1 

8+1 


Heq(s) “ ^2® 
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From (5) 


100 


K 


B + 20s + 100 


s + s + K(k^ + ^2®^ 


Equating the numerator and denominator polynomials, 

K = 100 

Kk^ » 100 

1+Kk2 >= 20 

Solving for and k^* 

K ° 100 

kj^ = 1 

k^ “ .19 

This design yields exactly the same transfer function as the 
Guilleraan-Truxal (GT) method. The two methods are not the same if a 
small change is made in the plant transfer ftcnction. The state variable 
feedback (SVF) system changes less than the GT system. 

If the sensitivity is defined as the percentage change in the 
closed loop transfer function for a small percentage change in a plant 
parameter, say K, then 

T' T'K^dKT 


where for simplicity 


T(s) 


Y(s)/R(s) 
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If K is the gain then 


T(s) 


KG(s) 

1 + KG(s) Heq(s) 


dT(s) G(s) K Gp(s) G(s) EeaCs) 

dK - 1 -h KG(s) He,(s) - ^ 


Using this result in the definition of sensitivity. 


S 


K 

T 


-^(gl 

KG(s) 


For the GT system, from (6) 


S 


K 

T 


8 (s+20) 

(s^ + 20s + 100) 


For the SVF system, from (6) 


( 6 ) 


S K „ sis+l) 

^ T 2 

a + 20s + 100 

For the GX system the sensitivity is greater at aLl frequencies than 
the SVF sensitivity especially at low frequencies. 

Modified Observer Systems 

[3 1 

Luenberger*^ ■* has shown that an estimate of the states in the 
plant can be implemented as shown in Fig. 3. When state variable- 
feedback is used from these states a specified closed loop response may 
be achieved, Muterspaugh^^^ has modified this configuration by block 
diagram manipulations to a form similar to that shown in Fig. 5. The 















Figure 4. Modified Observer Block Diagram 
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series coo^ensator shown in Fig. 4, C(s), reduces to the Guilleman-Truxal^^'^ 

compensation in the case where H(s) is one. Similarly this control system 

[ 2 ] 

reduces to the Heq type of control resulting from state variable 
feedback considerations when C(s) is one. Fig. 4 represents a more 
general configuration than either Guilleman-Truxal or Heq type compensa- 
tion. All three configurations sho\jn in the above figures are equivalent,, 
Assuming that only the output is available, the Heq configuration 
is not realizable since Heq has zeros but no poles . Guilleman-Truxal 
tjrpe compensation is sensitive to plant variations. Series and feedback 
compensation is used to achieve low sensitivity and realizability. If 
Heq has zeros in the LHP then the "modified observer" (MO) system, is 
slightly higher in sensitivity than Heq, but if Heq has zeros in the RHP, 
then the MO system is dramatically better. These comments are illus- 
trated in the case studies at the end of this report. 


Modified observer System Design 

This section derives the necessary conditions which C(s) and H(s) 
must satisfy such that the desired transfer function Y(s)/R(s) is 
achieved when the plant has a transfer function G(s). The relationship 
among these transfer functions is 


T(s) 


Y(a) 

xa M ■ 1 t f m 

R(s) 


' CCs) G(8) 

1 + C(s) G(s) H(s) 


(7) 


If the feedback signals shown in Fig. 3 are combined (using 
block diagram manipulations), the resulting transfer functions are 

A(8)=k2^j^Pj^(s)+k2 (f j^Pj^(s)p2(s)+/2P2<s))+. . .+k^(/^Pj^(s)+. . 


H(s)»k^+k2P^+k2p^p2+* . • *^N-1 
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Since the poles of H(s) and A(s) are arbitrary, the control systems shown 
in Fig. 5 have N-1 arbitrary parameters; however, the poles of A and H 
are identical. 

One more manipulation of the system shown in Fig. 3 results in 
the system of Fig. 4 where 


C(s) 


1 

1 + A(s) 


The zeros of C(s) are the same as the poles of A(s) ; therefore the N-1 
zeros of C(s) are also the same as the N-1 poles of H(s). 

From these facts (7) can be used to determine the denominator 
of C(s) and the numerator of H(s). Dropping the argument s in (7) and 
denoting the transfer functions by 

■ s 

- f 

- 1 

I(s) . I 


Eqn. 7 becomes 

CN GN 

Y ^ CD GD 

R “ CN 6N HN 

CD GD HD 


( 8 ) 


If the right hand is sinq)lified then 


Y CN GN HD 

R “ CD GD HD + CN GN HN 


(9) 


Since CN and HD are the same denote them both by F(s). Then (9) simpli- 


fies to 
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Y ^ GN P 

R “ CD GD + GN HN 


(10) 


cross multiplying on both sides results in the polynomial, 

Y(CD*GrH-GN*HN) - R*GN*P (11) 

Define the nmber of zeros in G as NG, the number of poles as 
ND, and similarly NY and NR for T(s), If the order of both sides of 
(11) is to be the same, then 


NY + ND - 1 + ND = NR + NG + ND - 1 

NY + HD •= NR + NG (12) 

NR - NY « ND - NG 


Thus the pole zero excess of the plant determines the pole zero excess 
of the desired closed loop transfer function. If like powers of (11) 
are equated, then NY + 2ND equations result. Since CD and HN are unde- 
termined the number of unknowns is 2ND. For NY « 0, the number of 

unknowns is equal to the number of equations. Just as in regular state 
f 2l 

variable feedback*' , for zeros in Y/R, additional states must be 
created. Let the order of P(s) by N-1, then (12) becomes 

NY + N- l + ND = NR + NG + N-1 

and for the nutrfjer equations to be equal to the unknowns, 

NY + ND + N 2N 


N »= NY + ND 


(13) 
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Thus, the order of both controller systems C(s) and H(s) must be one 
less than the nuni>er of poles of G plus the number of zeros of Y/R, 
assuming that (11) represents a set of linearly independent equations. 

The design procedure is as follows : 

1, Select the desired closed loop transfer function from the 
specifications . 

2, Select a polynomial, P(s). 

3, Equate the actual and desired transfer functions (10) and 
solve for the coefficients of CD(s) and HN(s). 

Program 1 

The first program (MODOBSl) confutes the transfer function of 
the series conq>ensator and the feedback compensator as a linear function 
of the denominator of the feedback compensator. The user must specify 
the plant transfer function and the desired closed loop transfer 
function. Thiis, this program specifies the design of the system except 
for the poles of the feedback compensator which the iiser is free to 
choose in any way he wishes. 

Program 2 

The second program (M0D0BS2) computes the transfer function of 
the series conq>ensator and the feedback compensator which minimize the 
low frequency sensitivity of the closed loop system to plant variations. 
The user must specify the plant transfer function and the desired closed 
loop transfer function. If the time constants in the plant are much 
lower than in the closed loop system, the poles of the feedback con^en- 
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sator may be in the right half plane (RHP). If the poles of the compen- 
sator are in the RHP then the system will be very sensitive to variations 
in the parameters of the compensators. In this case, it is better to 
use one of the other programs. This program is limited to systems with • 
no zeros. 

Program 3 

The third program (MODOBSS) computes the weighted integral 
sensitivity of the system for a given polynomial P(s). The user must 
specify the plant transfer function, closed loop transfer function, the 
weighting to be given the sensitivity as the frequency varies, and the 
denominator of the feedback transfer function, P(s). The weighting 
function must have more poles than zeros. 

Input Data 

In many circumstances the available data for the plant is in 
state variable form. This data is converted to a transfer function by 
the STVFDBK program. The details of the program are included in reference 
1. The program listing and input data format are also included in the 
Appendix. 

Output Data 

The program, B0DE4, is provided to plot the transfer functions 
conqmted in the above programs on a logarithmic scale. For each program 
some of the transfer functions of interest are: 

1. Plant transfer function 

2. Closed loop transfer function 
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, Sensitivity as a function of frequency 

4. Loop gain function 

The loop gain function divided by one plus the loop gain is the proportion 
of measurement noise at each frequency which will reach the output. Thus, 
if the loop gain is much greater than one then nearly ail the noise 
reaches the output, and if the loop gain is much less than one very- little 
measurement noise reaches the output. If the phase shift of the loop 
gain is nearly 180“ at a gain of one, the measurement noise will be 
amplified. 

Program 4 

The fourth program (M0D0BS4) is useful to check the results of 
the design. Given P(s), and the results of the first program (MODOBSl) 
this program con 5 )utes HS(s), CB(s) and the roots of both the series and 
feedback compensators. This program and the first program give a 
conq)lete determination of the design. 

Results 

A discussion of each program and several case studies are inclu- 
ded in the Appendices. MDDOBSl is useful to compute HN(s) and CD(s) in 
terms of the coefficients of P(s). Quite amazingly the coefficients of 
HN and CD are linearly dependent upon the coefficients of P. M0D0BS2 
confutes the coefficients of both compensators which achieve minimum low 
frequency sensitivity. Essentially this result is achieved by forcing 
CD(s) to be an N-1 order integrator which increases the gain at low 
frequency, and thereby decreases the sensitivity since = g (g) g (s ) • 
Unfortunately if the plant denominator coefficients are smaller than 
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those of the denominator of the closed loop transfer function the feedback 
compensators may have poles in the RHP. MODOBS3 determines the weighted 
integral sensitivity when P(s) is specified. M0D0BS4 is useful for 
checking the roots of the resulting compensators after a design is 
complete. Test results are given in the appendix for three examples. 

These examples show that an improvement in the system sensitivity over 
that of just a series compensator can be achieved even if only the output 
is available. When the closed loop response is much faster than the open 
loop, the sensitivity of the MODOBS system approaches the sensitivity of 
the Heq system (which is best but not realizable) if the roots of P(s) 
are made much larger than the closed loop roots. If the closed loop 
response is not much faster than the open loop response, then the MODOBS 
system may have a sensitivity lower than the Heq system. This situation 
exists whenever Heq(s) has zeros in the RHP. These results point out 
the important fact that a sacrifice in sensitivity to plant variations 
is made by increasing the system performance requirements. Because of 
this fact, careful attention to the plant performance should be made for 
lowest sensitivity and unnecessary demands on overall system performance 


should be avoided 
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Appendix I - MODOBSl 
Dls cussion 

The basic algorithm of this program is based upon (5). First the 
following polynomials in s are computed; 

T ■= (GD) (Y) 

Q = (GN) (R) 

V = (GN) (Y) 

so that (5) becomes 

(T)(CD) + (V)(HN) = (Q)(P) (Al) 

Each of the terms in (Al) are polynomials in s where P(s) is unspecified. 

2 

For instance, P(s) «= P^s + P^. The coefficients of CD(s) and HN(s) are 
the unknowns. 

As aii exan5>le, suppose the plant is second order, then P(s) is 
first order. In this case (Al) becomes, 

(CDaT^ + + (CD2T2 + CD^T^ + H2V2 + ^^^2)3^ + (CD 2 Tj^ + CDj^T 2 

+ H2V^ + H^V2)8 + (CD^T^ + H^V^) . P2Q3S2 + (p^Q^ + 

+ (P2Qj^ + PiQ2)s + P^Qj^ (A2) 

Equating coefficients of like powers of s results in the equations 

I3CD2 + V3H2 - Q3P2 

IjCD^ + I3CD3 + + V3H3 . QjP2 + Q3P3 
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T^CD^ + T2CD^ + = ^2 + Vl 

^l“l “ Ql^l (A3) 

These equations can be rewritten more coiqiactly in matrix notation as 



The matrix on the left will be defined as' matrix A, and the matrix on 
the right will be defined as matrix C. Then 



Therefore A relates the unknowns CD^ and to the parameters P^. 
MODOBSl forms the matrix A from the two arrays, T and V, and then inverts 
A. Next the C matrix is formed from the Q array and then A~^C is 
computed and printed. 



Ixxput Fonnat 


Card No. 


1 


Columns 


Description 


Fonnat 


1 

NY “ order of Y(s) in Y/R 

2 

NR = order of R(s) in Y/R 

3 

NG “ order of GN(s) in GN/GD 

4 ' 

ND = order of GD{s) In GN/GD 

5-80 

Identification of problem 

1-10 

Coefficient of in Y (s) 

11-20 

Coefficient of s^ in Y(s) 

etc. 


71-80 

Coefficient of in Y(s) 

1-10 

Coefficient of s^ in R(s) 

etc. 


1-10 

Coefficient of s^ in GN(s) 

etc. 


etc. 

-L 

Coefficient of in GD(s) 

j : — 


Notes : 

1. If more than eight coefficients are required then two cards may be 
used to identify that polynomial. 

2. A decimal point must be used for each coefficient, but not for card 
No. 1. The decimal point may be placed in any of the ten columns 
for each coefficient. 
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3. As many problems may be run as desired by repeating the above format 
for each problem. 

Output Format 

See exasp le. 

Subroutines to be tJ-sed 

PLYMLT Multiplies two polynomials together. 

MXINV Inverts matrix A. 

FSNXBAT Prints the ratio of two polynomials. 

Flow Diagram 













o 'n n N n n r) o n 
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PROGRAM MODOB51 (T MP UT ,0 U't PU T , T A P E '4- IK' PU T ) 

ODIMENSION YflOi R ( 1 0) r GD(10)» G-MflDt V(20)t T{20)» Q( 2G ) 
l.A(3Gi36)»Al(3Gf3P. )fC(36fl)J)»B(3 6tlfUiID'(6) 

Y(SJ/R(S) IS DESIRED CLOSED LOOP RESPONSE. 

GP(S)/GD(S) IS PLANT TRANSFER FUNCTION. 

NY IS CRDF.R OF NUH. IM Y(S)/R(S). NR IS ORDER OF DEN, 

NG IS order OF NUM. OF G(S). NO IS ORDER OF. DEN. . 

NR -NY MUST EQUAL ND-NG 
READ 101 »NY,NR .NG. NDtID 
PRINT 200. NY .NRt NG .ND» ID 
IF(E0F.4) 99.2 
NY Ir NY 1 

read 102. (Y (I ) . Ir 1, NY 1) 

Yfl) IS COEFF. of S*'.(I-1) in Y(S). 

NR 1-NR+l 

R{I) is COEFF., 0F_ S**< I-l). in. P{S) 

READ 102. IR(I)» Irl.NRl) 

CALL PRNTRATfY.NVl.R.NRl.lHT.lHS) 

GN(I) IS COEFF.. .OF^S’k^dr-l L.IN...GNtS)_., 

NGlrNE^l 
NOl-NO + 1 

.READ 102. fG.NCI). Irl.NGl) 

GO(I» IS COEFF. OF IN GDtS). 

5 read 102. (G0( I) .1 = 1 .NDl ) 

.CALL. PRNTRAT (GN, NG 1. GD .NOl .IHG tlHS.) 

C FORM CHARACTERISTIC POLYNOMIAL. 

CALL PLYMLT ( G D . NO 1 » Y . N Y I. T , 0) 

CALL- PLYMLT. . IGN. NGi. R. NR 1. Q, OJ 

CALL PLYMLT (GN . N G 1 . Y » N Y I . V .0 ) 

NYNG=NYi+NGl 

=NR+NG+i 

00 G IzNYNG.N 

C NA IS ORDER OF A MATRIX WHERE A IS MATRIX RELATING UNKNOWNS 

_6_ .V(IJ-0.. TO PAR AM,. P. 

NA= 2* {N-1 ) 

DO 20 Jrl.NA 

. IF_ ( J. G T. . (.N- lJJt.G 0_„T.0..„I.2 

DO 1 0 I rl . N A 
IND=N-1+J 

IF .{ IND....LT., l)_GO_TO_a 

IF (IND .GT.N) GO TO 8 
At I. J) = T( INO) 

SO. TO. .1D_ - 
.8 A ( I » J ) rO . 

10 CONTINUE 

G0._.T0.2D 

12 DO 20 I-l.NA 

IND =l-I+J 

IF. { IN D . LJ_.1,)„G0 ..J 0„i8— 

IF (INO. GT.N) GO TO 18 
At I. J) =V (INO) 

G0_T0'..2O.-. 

18 A(I.J>rO. 

.20 CONTINUE 

L- CALL MXINV. ( A.NA 

C C IS THE COEFF. ARR«r. 

NM IrN- 1 

_D0„24„I=.1.» MA 

DO 24 Jzl.NMl 
IN0-2*N-I- J 
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24 

C 

C 


30 


4 D' 


“so 

99 

"101 ■ 
102 
200 
'201 
202 
204 


3 . 


4 


e . 


7 


13 


14 


IF ( IND .L T . 1) GO 70 2 2 
IF ( IND.GT .N ) GO TO 22 
C( It J) - Q (IND) 

GO TO 24 
Cl It J) rO. 

CONT T\'UE 

B IS AI#-C. THE ROWS OF 5 aR E tOOflL TO CD AND HN . 
COLS. ARE COEFF, OF P PaRAPETERS. 

00 30 Irl.NA 

DO 30 K tNMl 
B { It K ) rO . 

DO 30JrltNA 

B(ItK)rAl{ItJ)*C(JtKJt-B(ItK) 

PRINT 20 It [K tKrl ,NH1 ) 

DO 40 IrltNHl 
NCrN-I 

PRINT 202tNC.{ 3( ItKl'tKzYtNMl )" 

DO 50 I-NtNA 
N3=NA*1-I 

'PR I N T 20 4 r' N B t ('6 ( i t K K = lt"NM 1 
GO TO 1 
STOP 

'FORM AT' ■( 4l lt6A"lO*)"’ 

FORMAT (8FID.3) 

FORMAT ( IHlt 4( ZX tl 1) t2Xt 6A 10 ) 

FORMAT ( IHOt 3X't8 (3'H' P( VI It" IH) t'OX V)" " 

FORMAT (4H C0( t I It 2H >= 1 1 0( El 0. 3t 3X ) ) 

FORMAT (4H HN C 1 1 1 1 2H )r , 1 0( El 0. 3. 3X ) ) 

END 

subroutine' PRNTR AT < A t NN tB tN Ot 0 EPV AR . INO V A R ) 
DIMENS ION AJ 20 ) f_B_( HJ 2 0 ) t N ( 2G ) t AT { 2 0) tBT ( 20 ) 

“ne":;!' . - * 

DO S IrltNN 
AT ( I )rO. 

'if { a' til > '4l6t 3 
H(NE ) = 1H + 

•AT (NE) =A (I) 

‘go TO 5 
H(NE)rlH- 

,_AT {NE1 -- A(IJ 

N{N£)zI-'i ' 

NErNEt i 
^CONTINUE 
NE’zNE-1 ■ 

PRINT 10 0 

” IN = M I_N 0 ( N E , 1 0. ), 

PR IN t“ i 0 i t ( N { I J Vi z 1 1 MIN)'" 

PRINT 10 2» (H (I I. AT (I J, INDV AR tizl tMlN ) 

IF (NE-10) 8t8.7 ' 

PRINT IDlt {NtlltirllTNEl 

PRINT' 10 2t (H(I )f AT (I ), INDVAR tin l.NE) 

„NF=i 

00 16 IzTtNO 
BT tl )zO. 

IF I B(I) 14»16ta.3„ 

H(NF )=1H + 

BT (NF> z3 (I ) 

GO TO ^15 

HtNF) zlH- 
BT(NF) 



15 N(NF)zI-l 
NF zNFi- I 

16 CONTINUE 
NFrNF- 1 
HHrlH- 

NHzMAXG( 13*NE» 13*NF) 

NHzMINO( 13Qf NH ) 

20 PRINT 103f DEPV AR ,INDVAR, (HH» Irl,NH)' 

MINzMINO (NFt 10 ) 

PRINT 10 1, (N (I ) , Iz 1, MIN) 

PRINT 10 2t (H a ) . BT (I )f IN Dv'fiR'f'l’zl tMIN )"’ 

IF (NF-10) 28»28t27 

27 PRINT 10It{N(I)» Izll.NF) 

PRINT 10 2. {HlIJ.BTm.TNDVAR .IzII.NfT 

28 CONTINUE 

100 FORMAT ( IHO) 

lOl' "format (6X'.lD( ilX.'l2 ))"“ 

1Q2 FORMAT ( 6X » 1 D ( A 1 . G 10 .4 .Al.lH )) 

103 format I IX »A 1. 1H ( . A1 »2H) 130A 1) 

RETURN 

END 

SUBROUTINE MXINV <R. N, Rl ) 

<- ‘"SUBROUTINE' TO F I NO T HE" I NV ER SE‘ 0 F ‘ a " G IYEN 
C MATRIX R. N is ORDER OF MATRIX. RI IS INVERSE 

^ matrix* subroutine uses GAUSS-JORDAN REDUCTI0N_» 

"C ■ R matrix 15 preserved. DIAGONAL ELEMENTS OF R'“MU5 

‘be nonzero'. 

_C jDIMENSION _R( 36 .36) »R A{ 36 .7 2) I RI(36.3^)_' _ 

C statements 20-26 ENTER R ARRAY INTO RA ARRAY 

■ C AND SET LAST N COLUMNS OF RA ARRAY TO IDENTITY 

C , MATRIX ... 

2 0 DO 2 6 I *r 1 , 'n 
21 DO 24 JzlvN 
-2 2__RA {_! .. J )_ _z R ( I , J) 

2 3 NJ z N + J 

24 RA(I.NJ) z D. 

25 NI z N + I 

26 RA(I.NI) z ’ 1 . 

C 

: ST ATEMENTS ..1_-1 2„RE.0UCE_ MATRIX. RA SO -THAT FIRST N 

C COLUMNS ARE SET EQUAL TO THE IDENTITY MATRIX 

1 NP z 2 * N 

2_D0. 12 I z 1 . N 

C - 

C STATEMENTS 3~5 ARE USED TO SET MAIN DIAGONAL 

_C ELEMENT J 0,. _U N I T Y_' 

"3 ALFA z RA{’f. i) 

4 00 5 J = I . NP 

5.RA(i,J) z RA(I.J) / ALFA 

C 

C statements 6-11 ARE USED TO SET ELEMENTS OF ITH 

_C „.COLUMN.„TO ..ZERO. 

6 00 1 1 K z 1 , N 

7 IF (K - I) 8. 11. 8 

8 beta r RAIK.I). 

9 DO 10 J : I. NP 

10 RA(K.J) z RA'(K.J) - beta ♦ RA(I»J) 

.1 l._CONTINUE„ 

12 CONTINUE 
C 



n o n o 


C 

C 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


STATEMENTS 3P-33 SET INVERSE MATRIX RI EQUAL TO- LAST 
N COLUMNS OF RA ARRAY 
30 DO 33 J-1*N 

3 1 JN - J + N 

32 DO 33 IrltN 

33 RKIfJ) r RA(I.JN) 

34 RETURN 
END 

SUBROUTINE PLYMLT ( A f L f Mt C. N) 

MULTIPLY ONE POL YNOMIAL BY"aNOTHER 

DEFINITION OF SYMBOLS I N AR GU MF NT L IS T 

A(I3» MULTIPLICAND COEFFICIENTS IN JHEOROER A(I) *S**(I-1) 
Lt NUMBER OF COEFFICIENTS OF A 

B(I)t MULTIPLIER COEFFICIENTS IN THE O RDE R 8 ( I) *S**_tIj-l 
Mt ■■ NUMBER OF COE FF I CIE NT S' CF 'B 

Cn)f PRODUCT COEFFICIENTS IN THE ORDER C (T? S* * ( I- 1) 

N. NUMBER 0F_ COEFFICIENT'^ '' 

REMARKS 

IF NzO» C(I) SET TO ZERO AND PRODUCT FORMED. OTHERVISE T HE 
[pr”o6uCT AND SUM NEWCr'OLD C '+ A*B IS FORMED. 


DIMENSION A(10)f B ( 1 0) ♦ __C (_20 )_ 

LPMrL^H-l’ 

IF (N) lOtlOflZ 
10_,__D0 11 A=1»LPM 
if ' C{'j'>-oVo 
12 DO 13 JrltLPH 

HAX-MAXO (J + l-M »1_) 

’MIN-MIND (L tj'V 

DO 13 I::MAXtMIN 
_1 3 _ _ C ( J ) zA ( I ) *B (J*1 ^I) ♦ C(J> 
"RETURN ' ’• 

END 


34 34 I n'co'RE' thermionic REACTOR 
1.16708 48.1,228 384. 520 200. 

.1.1670 8. 4 8 . 7 06 3 5 „4 u 8 . .5 3 15 ^39 2 .2.6. _ J_00 

1-16708 ' '48.1228 384. 520 200.-- 

5. D3585E-6 1. D542E- 3 . 0568639 1 ,8804 1 . 

END OF INFORMATION 
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Appendix 2 - M0D0BS2 


Discussion 

If CD(s) in (5) is defined by CD(s) = (Bl) 

where N is the order of the plants, then the system will have minimum 
low frequency sensitivity because the forward loop gain will be increased 
at low frequencies. If there are no zeros in Y(s)/R(s) or in G(s) . (5) 
reduces to a simple set of linear equations, shown below. For simpli- 
city and are set equal to one. For zero steady state error, 




C 


N-1 


^N-2 ^^N-] 


-N-2 


^1 ‘ ^2 ^ 3 ^ N -1 ^1 

• Rj^/Y 

HN^ =- . R^/Y + P^ 

= ^1 * + ^2 • + ... + 1 - GD^/Y 


Where the C terms are defined by 




S = S 


(B2) 


S-1 “ ■ ^} 
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The Eqns. (B2) are simply solved in the order shown fay hack 

substitution. That is, used. to solve for ^ which may he 

used to solve for P„ etc. 

N-3 


Input Format 


Card No. 

Co lumns 

Description 

Format 

1 

1 

N = Order of the Plant G(s) 

. 11 


2-80 

ID ■= Identification of problem 

8A10 

2 

1-10 

Y = Numerator of Y/R 

F10.3' 

•3 

1-10 

R(l) ~ coeff. of s^ in R(s) 

F10.3 


11-20 

R(2) = coeffc of s^ in R(s) 

F10.3 


etc. 


9 


71-80 

R(7) = coeff. of s^ in R(s) 

F10.3 

4 

1-10 

GN = nxamerator of G(s) 

F10.3 

5 

l-l'O 

• 

GD{1) = coeff. of s*^ in GD(s) 

1 « 

F10.3 


71-80 

* 

! GD(7) = coeff. of s^ in GD(s) 

F10.3 


See notes following input format for MODOBSl. 

Output Format 

See example. 

Subroutines to be Used 

PRNTRAT - Prints the ratio of two polynomials. 

PROOT - Determines the roots of P(s) to see if any roots are in 


RHP. 
















on* oj 1 rj M o o o n 


. . , - 

PROG^A^: M000P.S2 < I RP UT , 0 U'" PU T , T A PC 4 = TM T ) 

OirtRS lOM R ( 10 ) . GO ( 1 (i) fP { 1 CM t C ( 1 0) ,10(3) , HN ( 1 G ) , CO ( 1 3) •. 

C U ( 10 ) , V 1 1 0 ) 

C PROGRAM TO CALCULATE MINIMUM LOW-FREG SENSITIVITY, MODIFIED 

c observer system. 

C Y(S)/R(S) IS DESIRED CLOSED LOOP RESPONSE. 

C GN(S1/GD(S1 IS PLANT TRANSFER FUNCTION. 

C N IS ORDER OF 0£N. IN Y(S)/R(S), 

C N IS ORDER OF DEN. IN GN(S)/GD(S). 

C R(I) IS COEF. OF IN R(S). 

C 

C_ ASSUMPTIONS 

C *♦**#*+:►♦♦♦**♦ 

C TRANSFER FCN. GAINS ARE NORMALIZED SUCH THAT R(N + 1) = 1. 


C AND. GO (N+1 )zl. 

DESIRED TRANSFER FCN. HAS ZERO STEADY STATE ERROR. THUS 
i^ii*m******tt*<^*-*** Y : R ( 1 ) • 

GO(I) IS COEfT 6'f S* *( I-l) INGO (S ) . 

Y(S)=y » GN(S)rGN (NO ZEROS ALLOVEO). 

READ 101 . N,”i 6 
IF (EOF )99i Z 
READ IOZjX.. 

READ 1 D2 V'~{ R ( I ) , Iz 1 , N ) 

NP1=N+1 

GO ( N PI ) z 1. S R .{.NPIJ = !• 

READ 102 » GN 
READ 102 r (GDI I) tl = l .N ) 

PRINT 20D,Jl,iD 

'print' 201 ' 

CALL PRNTRAT ( Y t 1, R t NP 1 , 1H T. IHS ) 

PRINT 20 2 , 

CALL. PRNTRAT (GNtlVGD'iNPl '» 1 H6 »1HS) 

SOLVE FOR de n. OF H» P ( S). 

'nh‘izn-i - . — 

DO 40 I- If NHl 
IPlzI+1 __ 

‘ Tc ' C { I ) zGD ( IP 1 ) -R ( IPI ) 

P(N) zl . 

_00.5a Ir 
NI zN-I 
P( NI )zC( NI ) 

.IF(I.LE.l) G0..T0..5D-... 

DO 4 9 Jz2rl 
J2ZN-J + 2. 

,.JIrNI+.J-l 

4' P(NI )zP( NI )-R( J2 )*P( JI ) 

51 CONTINUE 

Call proot. ( nm.i, p, u,..Vf i)„.. 

PRINT 205 

PRINT 206» (U (I )f V( I) ► IzltNMl) 

-d 

C SOLVE FOR NUM OF H»HN(S) 

DO 6 0 'izlfN 
HNdlzD. . . 

DO 60 Jzl,I 
IJzI -J*l 



6C 


60 


33 .. 

101 

102 


200 

201 

202 

203 

204 

205 

ZQG 

210 


3 ' 


4 

■ 5 “ 

6’ 


7 

8 


13 


- 14 
. 15 


N,\n I )rP( J) (IJ) /Y *HN( I) .32 

HW (W )-HN (N )-GD a )/ Y 
PRINT 203 

CALL PRN TR AT (HNf N* P » Nf IHHt IHSJ 
DO 80 Irl. NMl 
CO ( I)zO. 

CD (N )- 1, . . . : 

PRINT 204 

CALL PRNTRAT (P ,N ,C0» N, 1HC» IHS) 

G r Y/GN . , 

PRINT 21 Ot G 
GO TO 1 

STOP 

format (I li 8A 1G ) 

FORMAT (3F10, 3) 


FORM AT tlHl tlOX tl If 4X .8A1 0) 

FORMAT (IHOf^THL DESIRED CLOSED LOOP TRANSFER FUNCTION IS*) 

FORM AT < 1 HO f*THE TRANSFER FUNCTION OF THE PLANT IS*) 

FORMAT tlHO f*TH£ TRANSFER F UNCT ION_OF , THE „FEFDB AC_K COMPENSAT OR IS*J 
FORM AT (IHO f*TH£"'TRANSFER FUNCTION OF THE SERIES COMPENSATOR IS*) 
FORMAT (iHOfSX* THE ROOTS OF P (S J ARt*13X*REAL P AR T* 10 X* IH A G I N AR Y 
PART*) 

format ( 32X»'2E20.’7 )' 

FORMAT UHO f*THE ADDITIONAL FORWARD LOOP GAIN REQUIRED IS 
1 Y/6Nr*Ein. 3) 

’ .END 

SUBROUTINE PRNTRAT ( At NN f 9 t N D f DE PV AR t IND VA R ) 

D I ME NS io N A (20 ) t B ( 20 ) »H { 20_) t N( 20)t AT (20 ) _»3J ( 2D L 

DO 6 1=1 fNN 

AT ( I ) = 0« „ _ _ _ _ 

IF' '('a (I) ') 4»’6t'3 
H(,NE )=:1H + 

AT(NE)=A(IJ 

60 to's"'"" 

H(NE ) = 1H- 
AT (NE) =- A ( I ) 

'N( NE'lrl- 1 
NE=NE<- 1 
CONTINUE 
NErNE- r 
PRINT 100 

_M I N = M IN □ { N E 1 LO. > 

’print lOlf {N(I)tI=lfMIN) 

PRINT 102» (H (I ) , AT tl )♦ INDV AR-tlfl »MlN ) 

IF (NE_-10) 3f8f7'„' ■ 

PRINT lOlf (N (I ) » 1 = 11 »Nb 

PRINT 10 2f (H (I ) » AT {I )f INOV AR .1=1 1. NE ) 

_NFrl_,_ 

DO 16 I=i»NU 
BT (I ) = 0« 

IF ( B(I) ). 14_fl6fl3„. 

~H( NF )-lH + 

8T (NF) z8 (I ) 

_G0 TO, 15. 

H(NF) =1H- 
BT (NF) =-B( I) 

„N(NF)rI-.l ... 

NF=NF*1 
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16 CONTINUE 

NF rNF - 1 
HHzlri- 

NHrN A.X D( 13*NE f I3*NF ) 

NH-MTNOC 13D»NH ) 

20 PRINT 103»DEPVAR»INDVAR. (HH» IrlfNH) 

MINzMINO (NF , 10 ) 

PRINT 10 1» (M (I ) . Izl , MIN) 

PRINT lG2i {H{I).3T(I)fINDVflR.Izl,MINl 
IF (NF-in) 28t28f27 

27 PRINT 101, (hid), inifNF) 

PRINT 102» (H(I)»ST(I),INOVARtIzil,NF] 

28 CONTINUE 

100. FORMAT (IHO) 

101 FORMAT ( OX ,ICH HX,I2)) 

102 FORMAT ( 6X » 1 D( A1 ,G 10 .4 »A1,1H )) 


103 FORMAT { IX t A 1 » IH ( » Al »2 H) 13 OA 1) 

RETURN 

END 

S U BROUTrNE.PROOT.(.N».Af.Ut R .) 

DIMENSION A( 20 ) »U( 20 )» V( 20 )♦ H( 21 )»8( 21 )t C( 21 ) 

IREVrIR 

NCzN + 1 

DOlIri .Nl, 

H( I) rA ( I ) 

„CONT.INUE_ 

PzO. 

0 = 0 . 

_ ,RzO. „ 

3 IF (H ( 1 ) ) 4, 2» 4 
2 NC=NC-1 

„V’( NC )-D. 

U(NC )-0. 

D01002I-1.NC 
H(I) =H(I+1 > 

1002 CONTINUE 
GO TO 3 

_4_IF (NC- 1) 5, lOD . 5 

5 IF (NC-2) 7,' 6f 7 

6 Rz-H(l)/H(2) 

GO TO 50. 

7 ' IF {NC-'3) Ot’St'g 

8 P=H( 2) /H (3 ) 

Q=H( 1) /H (3.)_ 

60 TO 70 

q IF ( ASS (H (NC- 1) /H (NO )- ABS( H( 2) /H ( 1 ) ) )1 Of 19 »1 9 

.LO_IREVz-IREV_ 

M=NC/2 

O0ilI=l,M 

NL-zNC-*- 1 =lJ 

F=H( NL ) 

H( NL )rH( I) 

.11 .Hi I) :F. - 

•IF(Q)13f 12fl3 
12 P-0. 

GOTO 15- 



13 PzP/0 
0= 1 . /O 

15 IF (R ) 1 G. 19 f 16 

16 R::l./R 
13 Er5*E-JC 

B < NC )-H ( NC ) 

C< NC l-HC NC ) 

3(NC+1 )r0, 

C( NC + 1')-G,. . , 

NP rNC- 1 

20 00 H 9 J- 1 » IGOC 
D 021 I 1 zl iNP. 

IzNC-Il 

B( I ) -H a ) + R*B( 1+ 1) 

..21..CtI) zB{I)*R*C(I+l) 

IF ( ABS (B (1 )/H( 1) )-E) 50 »5 0. 24 
24 IF (C (2 J ) 23 »22t 23 
-?■>„ RzR»J..„ 

GO TO 30 

23 R=R-8( 1) /C (2 ) 


30 D0.37Ilrl*NP 
IuNC-I 1 

B(I)rH{I)~P*B(I+l)-0 *6(1+2) 
,__X7 ..C ( I ) ( I .) - P * C ( I +_ 1 )„-0.*C.( I + 2 .)__ 

IF (H (2 ) ) 32 » 3 It 32 

31 IF (ABS (B (2)/H( 1) )-E) 33 »33» 34 

32 IF (ABS (B (2 )/H( 2) )- E ).33 » 3 3*. 34 

33 IF ( ABS (B (1 ) /H( 1) )-E) 70 .70» 34 

34 CBARZCI2 )-B( 2) 

D = C( 3) **2-C8AR;^C <41 
IF (0)36* 35*36 ' 

35 P=P-2. 

0r0*<0+l.>" 

G0T049 

36 PrP4-t8(2)*C<3)-8(l)*C(4) )/D 
0= Q + ( - a ( 2 ) ij: B A R+_B (_1 ) *C ( 3 n /D 

'49 CONTINUE ^ 

ErE* 10 . 

_G0T020 
^'50 NC zNC- 1 
V(NC)rO. 

IF ( IPE V) 51 .5 2* 52_ 

51 U( NC )zi. /R ■ 

G0T055 

52 U(NC)zR 

53 0054lr 1, Nc” 

H( I) zBtI+1 ) 

_54 CONTINUE. 

G0T04 

7 0 NCzNC-2 

IF (IREV) 71 *72*7Z. 

71 QPzl./Q 

PPrP/(Q*2. 0) 

G0T073 



PP -P /2.0 
73 Fr (PP) ♦♦ 2-OP 
IF (F jT'ii 75 »75 
7H U(NC + 1) = -PP 
. UCNC)r-PP 

V( NC + 1 )r SORT (-F) 

V( NC )r~V (NC + 1) 

G0T076 ... 

75 U(NC + 1 ) = -( PP/ABS (PP) )* ( ABS (PP) +5 0RT( F) ) 

Vt NC+1 )-0. 

U( NC )rGP /U (N.C+1) 

V( NC )rO. 

,76 D077TzltNC 
H(I)z8(I+2) 

77 CONTINUE 

■ G0T04 
LD.O..RETURN . 

END 

3. CASE STUDY II». THIRD ORDER._S.YS TEM_ 

1000 . . 

lODO. 200. 20. 

IDQD... 


2500 


100 
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Appendix 3 - M0DQBS3 


Discussion 

This program uses the following data: 

1, Plant transfer function G(s) 

2. Desired closed loop transfer function T(s) 

3. The denominator (P(s)) of the feedback compensator H(s) 

4, A desired weighting function W(s) on the integral over 
frequency of the sensitivity of the closed loop response 
to forward loop gain changes . 


The program then calculates the weighted integral sensitivity of 
the system. The method of calculation is based upon an algorithm by 
Effertz^^^ . The integrand of the integral sensitivity is calculated 
from 


2jtj . U 

W(s) 

S ^ (s) 

2 

ds 

-joo 

1 

-joo 

W(s)T(s) 


2rtj 

G(s)C(s) 


The algorithm is based upon an iterative manipulation of the coefficients 
in the numerator and denominator of the integrand in such a way that the 
integral of the squared magnitude results. 

A closed form solution for the integral 

1 = — i— r . . gp • _ f‘ F(p)F(-p) dP 

n 2rtj h(P)h(-P) . 2jtj 

(Bl) 

is presented in the paper by F. H. Effertz. The solution takes the form 
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of Equation 4 in his paper. A better algorithm for the computation of 
Equation 4 in Effertz is presented in the correspondence by Pazdera 
(Equation 7). A modified form of Pazdera 's algorithm has been coded in 
FORTRAN and is included. A test problem is, also, included. 

The use of Equation 4 in Effertz can be best illustrated by an 
example. Suppose we wish to evaluate the following 


It can be shown from residue theory that the correct answer is 


I 

n 


1 

s + c 


c - s 

2itj . 

-JOO 

2 

s +(a+b)s+ab 


2 . 

s -(a+b)sH-ab 


(B2) 


= - residue of F(s) F(-s) at s = a 
- residue of F(s) F(-s) at s = b 


(B3) 


where F(s) F(-s) indicates the function to be Integrated over s = Jm. 
Then the answer is 


2 2 
c -a 


2 2 
(c -b ) 


n 


2a(a+b)(a-b) 2b(a+b) (b~a) 


2 

c +ab 
2ab(a+b) 


(B4) 


To use Equation 4 from Effertz we must make the following associations, 
using the complex frequency variable s in place of p. 


g(P) 


= (s+c)(c-s) = 

/•r 


then 


P = s 
n-« 2 




2 

c 


2 

s 


h(P) 


s + (a+6) s + ab 


r 


P = s 
N = 2 


J a =1 
a^ = a + b 
a„ = ab 

V. 2 


(B5) 


then 



(B6) 


The integral^ for n — 2^ can be expressed as 


n 


Izil 

2a 


b b 
o 1 

a a 
o 2 


a. a 
1 o 

a a^ 
o 2 


2a 


Vz - °o'’ ll 

‘l‘2 


Substituting the correct values of a^, etc,, and cancelling minus 

signs for this problem we obtain 


T = ab -f 
n 2ab(a+b) 


(B7) 


What is most interesting is that the result depends only on the coeffi- 
cients of the known function being integrated and not on the poles of 
the function as one would suspect from residue theory. 

It has been mentioned that the algorithm from Pazdera (Equation 7) 
has been programmed. Some charges in the nomenclature of Equation 7 were 
necessary to facilitate coding. Note that the practice in Effertz and 
Pazdera is to make the highest membered coefficient correspond to the 
lot^est power of the variable p. Note, also, that subscripts such as a^ 
are in evidence. To facilitate coding, the lowest ,numbered coefficient 
was subscripted in the form A(l) and this coefficient was associated with 
the lowest power of the variable p, in this case p°. In general, the i*"^ 
coefficient A(I) is associated with the (i-1) power of p, P . (In 
FORTRAN the expression A(0) is not allowed.) Note from equation 1 above 
what must be done to operate on the function F(s) F(-s). If F(s) is of 
the form 

P(») = fg} (B8) 


then the algorithms suggested in both Effertz and Pazdera require the use 



of 
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g(s) = C(S) C(-S) 
Ms) = A<S) 


<B9) 


It is more desirable to input C(S) and A(S) rather than C(S) C(-S) and 
A(S). Consequently one uses the INTSQ subroutine by coding the coeffi- 
cients of C(S) and A(S) with ascending subscripts corresponding to 
ascending powers of S. The subroutine provides the operation C(S) G(-S), 
Furthermore, the subroutine checks to see If the lowest coefficient of 
C(S) or A(S) is zero so that factors of S may be either considered or 
cancelled. Two basic requirements must be met by A(S) and C{S). First, 
the roots of A(S) and C(S) must at least have non-negative real parts. 
(Hurwitz polynomial requirement) Secondly, the highest power of C(S) 
must be at least one less than the highest power of A(S) for convergence 
to be assured. If A(S) is Kth order, one inputs the S coefficients of 
A(S) and the N-l coefficients of C(S), 

The use of this subroutine INTSQ may be best illustrated by an 
example. Suppose we desired to evaluate the integral 

4 - s 
3 2 

4s + 3s - 2s + 

(BIO) 

The polynomials C(S) and A(S) as defined in Equation BS are input into 
the program as 

A(l) = 1 
A(2) 2 

A(3) = 3 
A(4) =* 4 



= J 


s 4- 4 

4s^ + 3s^ + 2s 4- 1 


C(l) = 4 
C(2) « 1 
C(3) « 0 


(BID 
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The program then utilizes the modified Pazdera algorithm and prints out 
the message 

THE VALUE OF THE INTEGRAL IS 12.250 


To check this answer, one can use Effertz Equation 4 for a third order 
case (n=3). 


= 


n 


2a 


•■o >=1 *>21 
% ^ ° 


ai aj 0 

% “2 ° 
0 


1 ”o“2V'2VrV3'’l 


2a 




(B12) 


For our example problem, the following associations follow 


b = 0 
o 



= 16 


Substituting the values of Equat 
computer's solution => 12.25. 


a = 4 
o 

^1 ^ (B13) 



B13 into Equation B12 verifies the 


Input Format 


Card No. 

Columns 

Description 

Format ^ 

• 

• 

See forma 

,t for MODOBSl 


« 

K 

1-10 

P(l) - coeff, of s° in P(s) 

8E10.3 


11-20 

P(2) - coeff. of s' in P(s) 

8E10.3 


etc. 





42 


Input Eormat (cont'd) 


Card No. 

1 

Coltunns 

Description 

Format 

K+l 

1 

NW - Number of coeff. in denom. of W(s) 

11 

K+2 

I-IO 

l'JN(l) - coeff, of in num. of W(s) 

8E10.3 



I’JN(2) - coeff. of s' in num. of W(s) 

8E10.3 

K+3 

1-10 

WD(1) - coeff. of s° in den. of W(s) 

8E10.3 


11-20 

etc. 

WD(2) - coeff. of s' in den. of W(s) 

8E10.3 

1 

! 



c 

1 

c 


.c 

c 


_a. 


10 


102 

103 

_Z01 

202 

20M 

.205 


C 

C 

C 

c. 

c 

c 
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PROGRAM HODOGS 3 ( INPUT, OUTPUT ,TAPE4r INPUT ) 

DIMENSION Yt lC)tR(10)fGD{10) tGN(10)iP{lD).CDflO) »HN(i0)fSN(40) t 
1 SO (4 0) »B (36t 18 ) » WN ( 1 0) t W 0 1 10 ) . SN T ( 40 ) t SO T ( 40 ) 

CALL MOOOBSl ( Y »R »GN, GO »5 »N Y1 »NR1 ,NG 1 »ND1 J 

N-1 IS THE NUMBER OF COEF, IN P { S ) »HN ( S ) . C 0 ( S ) . 

NrNR 1 + NG 1- 1 
Nl rN -1 

READ 1 02 1 1 P ( I) , Izl ,N II 

V(S) r WN/H'O must HAVE MORE POLES THAN ZEROS ... . „ . 

NW IS THE number OF TERMS IN THE DENOMINATOR OF W(S). 

READ 103.NW 

.READ. 102. (WN( 1) ,1=1 .NW)_ 

READ 102 , (worn ,lrl ,NW) 

PRINT 2C5 

CALL .PRNTRAT { V(N» NW ,W D..NW ,1 HW ..l.HS.) 

CALL PLYMLT(GDtNDl »Y »NY1 tSNf 0) 

N2=ND1+NY1-1 

_C A L L_ . P L Y K L.T ( VN N W .,_S N,t N 2 t„SN.T *, D > 

CALL PLYMLT(GN,N6I .R tNRI .SDt 0) 

CALL PLYMLTtWOtNW.SDtNZrSOT.C) 

N3zNV+N2-l 

N4 rN3+Nl ~1 
N2r2*Nl+l 
„D0 ID Jrl.JNl 
CD (I JrD. 

HN ( I )rD. 

.DO ID J-1»N1 

CD(I) = B(N-IfJ)=*P(J) + CO{I) 

HN(I )r3(N2~I»J)*P( J) +HN(I) 

_PRJNJ._2D 1 

CALL PRNTRAT (P»Nl»CDtNl, IHCt IHS) 

PRINT 2D2 

.CALL..PRNTRAT.(HNtNl ,P .N1,.1HH» 1H,SI 
CALL PLYHLT((:D,N1,5NT,N3,5N,C) 

CALL PLYMLTf ■ P,N1, SOT, N3 fSOf 0) . 

-N5rN4 ---1 

CALL INTSQ (S0.5N»N4» 5) 

PRINT 20 4» S 

L.GO. T-O.J. 

FORMAT (8 E10« 3) 

FORMAT III) 

„F0 RM AT.L*.Q_THE_TR.AN5F-£R— fLUN CT-IO-N_OFL-J-HE_S £iL.r£-S—COWB.EAl.SA-IOJ?_I3-t) 

FORMATING THE TRANSFER FUNCTION OF THE FEEDBACK COMPENSATOR IS*) 
FORMATING THE WEIGHTED INTEGRAL SENSITIVITY IS*E1D.3) 

_FORM.AJ„( *D„TJH£- WEXGHXLNG- F UNCX lOhL-OF— THE - IMTEGRAL— S£NS.I1L1-V IXY— ISj 
1 ) 

END 

..SUBROUTINE! .HODOBSl LY.tR ,GN, 60 ,B »N.Yl-fNRl-,N6l-.N04X 
□ DIMENSION Y(1Q)» RUG). GDUO). GN(IC). V(20). T(20). QC20) 
l.A(36. 36 ).Al(36. 36)tC( 3S.18) .8(36»l3)fl0(6) 

„Y(5)/R(S) IS DESIRED 'CLOSED LOOP RESTuu^f.. - 
GPIS)/G0(S) IS PLANT TRANS FER FUNCTION . 

NY IS ORDER OF NUM. IN Y(S)/R(S). NR IS ORDER OF DEN. 

_NG._I S OR DE.R.„OF„N Utl..._OF_ GXSJ._RD-US-0JlD£EL_QE_DEi4 
NR-NY MUST EQUAL NO-NG 
Yd) IS COEFF. OF 5**(I-1> IN Y(S), 
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f?(n IS COEFF. OF '?(Si. 

6N(I) IS COEFF. OF S**(I-1) IN GN’(S), 

GD(I) IS COEFF. OF 5**{I-1) IN GOCS). 

READ lUl »NY»NR .NGrNO.IO 
IF(E0F,4) S9»2 
NY IrNY +1 
NRl-NR+1 

READ 102» (Y II ) » I=lf NY 1) 

READ 102. ( R( 1) . I r 1. NR 1) 

PRINT 200. NY .NR. NG .ND. ID 
CALL PRNTRAT(Y,NY1»R.NR1»1HT.1HS) 

NGlzNG +1 
NOlzNO+1 

READ 102 . (GNl I) . Irl. NG 1) 

READ 1 D2 . ( GDI I ) .1 = 1 .NDl ) 

CALL PRNTRAT(GN.NGl.GD.NDl.lHG.lHS) 

FORM characteristic POLYNOMIAL. 

CALL PLYMLT IGD. ND i. Y. NY It T. 0) 

CALL PLYMLT ( GN » NG i . R . NR 1. 0. D) 

CALL PLYMLT IGN .N G1 . Y tN Y1 tV'.O ) 

NYNGrNYl+NGl 
N =NR+NG + 1 

"DO G IrNYNG.N 

NA IS ORDER OF A MATRIX WHERE A IS MATRIX RELATING UNKNOWNS TO 
V(I)=0. PARAM.t P. 

NA= 2* CN-1 ) 

DO 20 Jrl.NA 

.IF ( J. GT . ( N- 1 ) ) G 0_ JO„.l 2, 

DO 10 Irl.NA 
INDzN- 1 + J 

IF. ( IND -LT.l) G0_T0_8... 

IF (IND .6T.N) GO TO 8 
A{ I. J) ~ Tl IND) 

GO TO 10__ 

' A( I.’j) =0V 
CONTINUE 
GO TO 20 
DO 2D I=1.NA 
IND =1-I+J 

. IF { IND.LT.l ).. GO. TO. .18. 

IF (IND.GTtN) GO TO 18 
All. J) rV (IND) 

..GO TO 20 
A( I. J) =0 . 

CONTINUE 

call MXINV_ ( AtNA_»AI_). 

*7 'C IS THE COEFF. ARRAY. 

NH1=N- 1 
DO 24 IrltNA 
'do 24 Jrl.NMl 
IND=2>»N-I-J 

..IF (INO.LT. 1) GO TO 22_ 

IF (IND.GT.N) 60 TO 22 
C { I. J) = Q lIND ) 

GO TO 24 
C ( 1 1 J ) =0 . 

CONT INUE 

B IS Al-»C.THE Rows UF B AH-E ’EQUAL lU CU ANU HN . T H E CO LS_ .__A_R.E.. 
COEFF. OF P. PARAME TERS , 

DO 30 Irl.NA 



30 


40 


50 

..99 

IDl 
102 
.200 
201 
202 
2 0.4 



_C_ 

C 

C 

,.C. 

C 

c 

_c.„ 

c 

c 


_1C 

1 ] 

u 


. 1 .: 


_3. 


4. 

5 


DO 30 K zl tNMl 
B ( If K ) =0 . 

DO 30JrltNA 

PR INT 20 1» (K flCrl .NMl ) 

DO 40 in. NMl 
NC rN -I 

PRINT 20 2. NC . { B( I »K ) tKrl .NMl ) 

DO 50 IrN.NA 
N6zN&+ 1- I 

PRINT 204. NB. (8 ( I .K J , Kr 1. NM 1 ) 

RETURN 

9T0P 

FORMAT {4I1.6A1G) 

FORMAT {BF10.3) 

FORMAT ( IHl. 4( 2X .11) .2X. 6A 10.) 

FORMAT (10X.3(3H P C . 1 1 . 1 H) .8 X ) ) 

FORMAT (4H CD( . I 1. 2H )- . 1 0( El 0. 3. 3X ) J 
FORMAT (4H HN ( . I 1. 2H )= » 1 0( El 0..3..3X) L 
END 

SUBROUTINE' PLYMLT ( A. L. B» H. C » N ) 

MULYiPLY’~0 N e’ P 0LYN0VlMl”'BV'XN0’fHER 

DEFINITION OF SYMBOLS IN ARGUMENT LIST._ 

At I) MULTIPLICAND COEFFICIENTS IN THEOROER A(I) =»’s**(I-r> 

L. NUMBER OF COEFFICIENTS OF A 

» MULTIPLIER XO E£F_I.CI£N.TS.,.lN„IH£L_.0RJ3EE_BJ.JJ_*5jMl_lJLrJJI 

M. number of COEFFICIENTS OF B 

.CCI)» product COEFFICIENTS IN THE ORDER C ( I ) * S* * ( I- 1) 


remarks 

__IF„fcL:0 t.„.C( JJ_S FT T O 7FR 0 AND PRODU CT FO RMED^ OlilEJR WISE.. .IH.E_ER-QJ1IJC 

AND SUM NEWCz OLD C + A*9 IS FORMED. 

. DIME NS ION A( 10 ) »_,.B.(IDIi__C.L20J 
LPM=L+M-1 
IF (N) 10.10.12 

_D0. II „ J,=l .LPM 

C(J)rO.O 
DO 13 J=1 .LPM 
..MAXrMAXO(J+l-M.._L.) 

MINrMINO (L »J > 

DO 13 IrMAX.MI 
,C( J).rA a 
RETURN 
END ' 

..SUBROUTINE, P.RNTR.AT_ ( A. NN .B rND ..DEP.V AR.»1DD.V,A RJ 

DIMENSION A(20).B(20).H(20 ).N( 20).AT(20) .8T(2D ) 

NE:rl 

_D0 .6..I_rl.»NN 

AT ( I )zO. 

IF C A(I ) > 4.6. 3 

_H( NE )= lH-». 

AT tNE) =A (I ) 

GO TO 5‘ 

_H (.NE.)z,lH*: 

AT tNE) =- At I) 

NtNE )rl- 1 
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N£ 'NE+ 1 • 

6 CONTINUE 
NErNE- 1 
PRINT ion 
HINzMINC (NEf 10 ) 

PRINT ID 1 , f N Cl ) , T- 1 . MIN) 

PRINT lOZt CH (I ) t AT < I ) ♦ INDV AR »I =1 »M IN ) 

IF CNE-in) 3»3t7 

7 PRINT lOlf CN (I ) ♦ Irll .NE) 

PR INT 10 2» (H (I ) , AT ( I ) , INDV AR f T zl 1. NE ) 

8 NF=l 

DO 16 IrlfND 
BT (I )r:0, 

IF ( 8CI) ) 14tlcVl3 

13 H(NF)rlH + 

8T CNF) zB ( I ) 

GO TO 15 

14 HCNF) rlH- 

__ BT(NF)=-B(i: 

15 NCNFJ-I-'l 
NFzNF+ 1 

16 CONTINUE 
NF -NF- l“" 

HHzlH- 

NHzMAXCC 13*NEf 13*NF» 

NH-MTNOl 13D* NH ) 

20 PRINT lD3f DEPV AR tINOVARt CHH» IrltNH) 

HINrMINO (NFf 10 ) 

“print 1 0 1 » ( N ( I ~)"»"l z 1 »' M I N )' 

PRINT 10 2* <H(I ),BT(I). INDVAR fizl »MIN) 

IF (NF-IO) 23.28»27 

'27 PRINT 10 1. (N (I), inVtNF) 

PR INT 10 2, (H (I ) , BT tl ), INDV AR ,I zl 1,NF ) 

28 CONTINUE 

l6d“ FORM AfTlHoT" 

101 FORMAT { 6X ,10( 11X.I2 )) 

_102 FORMAT ( 6X , 1 0( A1 ,G 10 ,4 , A1 ,1 H )) 

103' format C IX »A1,‘ih(, AlV2H)‘zfY3 0Air‘ 

RETURN 

END 

SUBROUTINE MXINV (R, N ». ,RI ) 

C SUBROUTINE TO FIND THE INVERSE OF A GIVEN 

C MATRIX R, N' IS ORDER OF, MATRIX, RI ,IS INVERSE... 

'c' matrix’.’ subroutine’ USES GAUSS- JORDAN REDUCTION, 

C R MATRIX IS PRESERVED. OLAGONAL ELEMENTS OF R MUST BE NONZERO. 

D I ME NS 10 N_ R t 36 » 3 6 ) ,RA ( 36 ,7 2) , _RJJ. 3&*36) 

c ' 

C statements 20-26 ENTER R ARRAY INTO RA ARRAY 

C AND SET LAST N COLUMNS OF RA ARRAY T O IDENTITY 

"c matrix' 

20 00 26 Izl.N 
_ 21. DO 2 4 J z 1 
22 RA (I ,J) z R(I, J) 

2 3 NJ z N J 

24 RA ( I ,N J) ,z 0 . .. .. 

25 NI = N ♦ I 

26 RA ( I ,N I ) z 1 . 

_C . .... 

C" ” St'‘aTEMe'N*TS” 1 -IZ REDUCE MAT'rI)! RA SO THAT FIRST N 

C COLUMNS ARE SET EQUAL TO THE IDENTITY MATRIX 
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1 NP n 2 * N 

2‘ DO 1 2 I - 1 f N 

C STATEMENTS 3-5 Af?E USED TO SET MAIN DIAGONAL 

C ELEMENT TO UNITY 

3 ALFA r R A< It IJ 

4 DO 5 J r It NP 

5 RA(ItJ) r RAdtJI / ALFA 

STATEMENTS 6-11. ARE USED . TO , SET ELEMENTS OF ITM 
COLUMN TO ZERO 
60011 K-ltN 

7 IF (K-I) Stilt 8 

8 BETA = RAlKf I) 

3 DO 10 J = It NP 

10 RATKtJ) r RA(KtJ) r..BETA_* R A ( It.J) 

11 CONTINUE 

12 CONTINUE . 

STATEMENTS 30- 33 SET I NV ER SE^m'^TRiV RI~ E QU AL~ T O^LaTt" 

N COLUMNS OF RA ARRAY 

3 □ DO 3 3 J = 1_» _N 

Sl’ jN 'r j”+ N 

32 00 33 I=1»N 

35 RI (I tJ>. r RA (I_tJN)_ 

3 4 RETURN ' ' 

END 

SUBROUTINE PLY SO, (Ct.Nt BJ 

DIMENSION C(40)tB(4D) 

RETURNS B(S) - C(S)*C(-S) TO MAIN PROGRAM, 

B(I) IS. COEF.. OF..S*^Irl IN 3.(S). 

ClI) IS COEF. OF S**I-1 IN C(S). 

N-1 IS number of COEF, IN B(S) AND CIS), 

N1 =(.N-ll./2. 

00 20 IrltNl 
MO =-l 
_ B { I } :0 . 

II r 2*1-1 
DO 20 K=lt II 

MO r -1*M0 

20 B( D iB (I )*MO*C (K )*C( 2*I-k> 

NMlzN-1 

N2=N1+1 , . . 

DO 30 I=N2tNMl 
IINz2* I-NHl 

B ( I ) rO - 

M0 = (-1) ** NMl 
DO 30 KrIINtNMl 

HO = . „„ 

30 B( I) rB (I )*HO*C (K )*C( 2* I-K) 

RETURN 
END .. . 

SUBROUTINE INT SO ( A t C tN ,S ). 

DIMENSION B( 10 ) t At 11 )t C( 10 ) 

_C RETURNS SzINTEGRAL OF. C ( S) *C t- SJ /A tS »*.A{ -rS K._T.O H AI.N_EROGRAM,.. 

C BtS) HAS N-1 TERMS. AtS) HAS N TERMS. 

C BID IS COEF. OF S*»(2I-2) 

..C.^ Ctl) rs'.COEF. OF S**(It1J. 

C IF THE LOWER ORDER DEN. AND NUM. COEFF, ARE SHALL (LESS THAN 

C 0) THEN DIVIDE BOTH NUM. AND DEN, BY S**2. 
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0-1. OE -G 

1 IF { &8S{ ftt 1) >-0) 2 »2 r2D 

2 IF ( AISC Ct 1) )-D) 4 .4 .20 

4 IF CK-1) 5.5.6 

5 t^RINT 200 

6 PRINT 20 1. A{ 1) . C( 1) 

KrK V 1 

NrN- 1 

00 10 1:1. N 

A( I) zft (i+n 
IF(I-N) 8.10.10 
.8 C( I ) =C ( I +1 ) 

10 CONTINUE 

60 TO 1 

20 CALL PLYSO I C.N. 8) 

NM2:N-2 

DO 50 Krl.NM2 

NK:N-K 

■ BA:8(NK) /A(NK+1)’ 

.AArO . 

IF CK .EO. 1) GO TO 4_0_ 

‘aA:A (NK + 2) /A iNK+ 1) 

40 NMK:{N-K)/2 

DO 50 1:1. NMK 
NKI:NK-2*I + 1 

B(NK-I):B(NK-I )-BA*A(NKl) 

. A ( NK 1+ 1) :A { NKI +1.)- AA,»A_{ NKIJ 

50 CONTINUE 

S-B< H /( 2.*At2)*A{ 1) ) 

...... return .. .... 

200 FORMAT ( *0 THE FOLLOWI.NG .COEFF. WERE FOUND TO BE SHALL AND CANCEH 
lED IN THE NUMERATOR AND DENOMINATOR*) 

20 1 FO RM AT_ J_* T HE JDEN_P M_1N A TOR C O EFF.:* EID.3 * THE NU HERA TOR. COEFF. 

1:* E 10.3) 

END 

0303 ~C'aSE STUDY”lY“fHlW’ORDER~S>YfEH 
1000. 

..10 DO . 20 0. 20 l._ 

1000. 

2500. 100. 1, 

. 80 .- 1 ,^ 


.70.0 

? 
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Appendix 4 - M0D0BS4 
Discussion 

M0D0BS4 uses the matrix relating the coefficients of CD(s) and 
HN(s) to the coefficients of P(s) and also the coefficients of P(s) to 
determine both the series compensators and the feedback compensators. 
The roots of the compensators are included. 

Input foimiat 


Card No. 

Coltanns 

Description 

Format 

1 

1 

N » order of plant G(s) 

11 


2-80 

ID = Identification 

8A10 

2 

1-10 

CD(1, 1) coeff of array from MODOBSl 

E10.3 


11-20 

CD (1,2) coeff of array from MODOBSl 

EIO.3 

3 

1-10 

CD(2,1) coeff of array from MODOBSl 

EIO.3 

• 

• 

• 

11-20 

CD (2, 2) coeff of array from MODOBSl 

EIO.3 

N+2 

1-10 

HN(1, 1) coeff of array from MODOBSl 

EIO.3 

• 

• 

• 

11-20 

HN(1,2) coeff of array from MODOBSl 

EIO.3 

2N+2 


P(l) - coeff of s^ in P(s) 

EIO.3 


m 

P(2) - coeff of s' in P(s) 

EIO.3 






















FUNS S > 
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PROGRAM M00005 4 ( INPUT F OUTPUT ,T APF.4Z INPUT ) 

C PROGRAM TO CALCULATE C 0*^ PE NS A T OR S WITH P(S) AND RESULTS OF HOOOBSl 

DIMENSION ID S8 ) F CD ( 1 Rf 18 ) r HN a 8* IS ) r C( 18 ) » H 1 13 Jf U{ 23 ) * VS 2Q) ♦ P.t2Q ) - 
C N IS THE order OF THE PLANT, GSS) AND NO, OF COEF IN P(S), 

1 . READ 1 DO * N , ID 

IF(E0F,4) 99f2 

2 PRINT 20n» N. ID 
DO I n I -1 ,N 

ID READ ini. (COS I» J) 1» N) 

DO 15 Irl,N 

15 READ lOlf (HNUtJ), Jrl»N) 


read -1 01 . (P (I ) F Izl fN ) 

DO 20 I=1»N 

C ( I ) rO . 

H ( I ) rO , 

DO 2D JrlfN 

C( I) rCd l + CO (I fJ )*P{ J) 

HSI) rH ( I J + HN (I fJJ*P(.J). . 

PRINT 20 4 

CALL PRNTR AT (P fN fC fN f1 HC f1 H5 > 

_PRINT .20 5 

CALL PRNTRAT ( H f N f P , N» IH Ht IH S ) 

NH1=N-1 

.CALL PROOT. J NH.lF.C.F.U».V.t..I).., 

PRINT 201 

PRINT 202* (U(I J F V( I) .Izl »NMX J 

.CALL P..R 0 O.TJ.N M l_F.P.t U_f..V_f 1_) 

PRINT 203 

PRINT 202* (U (I ) f vm Firl fNHI J 

CALL PROOT INMltHN_F.y.F.V.FlJ 

PRINT 220 

PRINT 202f (Utl »* Vf I) fI:1 »NM1 ) 

GO.. T_0__1._ 

yy STOP 

100 FORMAT (IltSAlD} 

„10.1.. _FORMAT_ ( 8E10.3) 

200 FORMAT ( IH 1 , 2X f 1 1 f 2X *8 A1 OJ 

201 .format (1H0f5X*THE ROOTS OF C(S) ARE*13X*REAL PA RT * I OX >* I M AG IN A RY 

.PART.*} ^ 

202 FORMAT (32Xf2E20.7) 

203 FORMAT (1H0f 5X*THE ROOTS OF P(S) A RE *•! 3X *R E A L PAR T * 1 OX n«u IN A RY 

220 FORMAT (1HQf5X*THE ROOTS OF- HN (S ) ARE*13X*REAL P AR T* 10 X* IM A G I N AR Y 

•IPART*} 

_20.4 F 0 RM AT.„{_1H 0 f_* T.HE—TJ? ANS FE-R...FU N C.TXO N., 0,E_J.H £.^Ei?.I.ES— C.0 'M pE-NSAT DR_ g S J lX, 

205 . FORMAT (IHDr *THE TRANSFER FUNCTION OF THE FEEDBACK COMPENSATOR IS 
1*) 

.END 

SUBROUTINE PRNTRAT ( A f NN fB fN D f OE PV AR » IND VA R) 

01 ME NS ION A ( 20 ) * B( 20 )f H( 20 ) F N( J . AT .nTf 7n » 

-NE.rl.. 

DO G IzIfNN 
AT ( I)::D. 

IF A ( I } . 1. U. 
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3 H(NE)rlH+ 

AT (NE» rA (I ) 

GO TO 5 

A H( NE TH- 

AT (NE) r- A( I ) 

5 N(NE)zI-l 
HE =NE+ 1 

6 CONTINUE 
NE ~NE- 1 
PRINT 100 
MIN-MINQ (NE » 10 ) 

PRINT 10 1 . ( N (I ) , Ir 1, MIN) 

PRINT lD2t (H (I ) ♦ AT (I ) , INDV AR tI tM IN ) 
IF (NE-in) 8taf7 

7.. -.PRINT 10 1» (N (I-J » I- 11 fNE) - - 

PRINT 10 2f CH a ) , AT ( I ), INDV AR ,I rl If NE ) 
8 NFri 


DO 16 I-1»ND 
BT (I >zO. 

IF { B(I) J m»lGfl3 
_13„ ,.H(NF)zlH + 

BT (NF) z9 (I ) ' 

60 TO 15 

.14 _ H(NF) zlH- 

8T (NF) =-B( 1> 

15 N(NF)rI-l 

_NF-NF+1 

Tg”’ continue" 

NFzNF-.i 
HHrlH- 

NH=MAXO( 13*NE» 13-»NF) 

NHrHINO( 130fNH) 

,.?Q_ „PR IN T 10 3f DE PV AR.f I NOVA R« (HH» Irl ,NH) 

MIN-MINO (NF f 10 ) 

PRINT IDlf (Nd )f Irl, MIN) 

,,.PR I N T 1 0.2 f .( H ( I ) f BT ( I ) , 1 N.D^1R_, I r 1_, M I 
IF (NF-in) 28 *28, 27 

27 PRINT 101t(N(I), inifNF) 

_^PRINT 10 2, CH (I ),_8T (I^>, INDV AR ,Ir l 1,NF ) 

28 CONTINUE 

10[ FORMAT (IHD) 

101 FORMAT ( 6X ,.10( IIX, 12 )) 

102 ‘ FORiMAT ( 6X , 10( A1 ,G 10 ,4 ,A1,1 h ’)T 

103 FORMAT ( IX , A 1 , IH ( , A I ,2 H ) r, 13 OA 1 ) 

„RETuRN. _ 

END 

SUBROUTINEPROOK N, At U, V, IR ) 

DIMENSION A( 20).t U( 20.) ♦ V ( 20J.t..H.{.2U t b 7(.21 ).» C.( 21 
IREVrIR 
NCrN + 1 
„ D01Irl.,NC_ 

H( I) =A (I ) 



1 continue: 

p-o. 

Or 0 . 

Rr 0. 

3 IF ( H ( 1 n 4» 2f *4 
2 N'CzNC-1 

V ( NC ) = 0. 

. UtNOrO. 

DO icrzl: It NC 
H( I) zH (I+l ) 

.1002 CONTINUE . . „L 
GOTO 3 

4 IF (NC- U 5. lOOt 5 

5 IF (NC- 2) 7t G» 7 

6 Rz -H ( 1 )/H( 2) 

60T050 

.„7 IF (NC-3) 9t 8i.9.... 
e PzH( 2) /H (3 ► 

QzHU)/H(3) 

; _ GOT070 ^ 1 : 

9 IF(ABS(H(NC-1)/H(NC) »-ABS(H(2)/H(lH )10tl9*19 

10 IREVz-IREV 

_MzNC/Z 

DOllIzltM 
NL zNC^-l-I 


FzHI NL > 

H(NL )zH( I) 

11 H( I ) zF 

IF.(Q )13»„12jL3. 

12 P=D< 

GOTO 15 

...13 PzP/0 . 

Ozl. /O 

15 IF (R )16t 19.16 
_1.6„Rz: 1 , /R 

19 EzS.E-lO 

B ( NC )zH( NC ) 

C( NC )zH( NC )_ 

B( NC+1 )zO. 

C{ NC + 1 )zO. 

NPzNC-1, 

20 DOqgjz 1, 1000 
00 2.1Ilz:l ,MP 

_ ... IzNC-I.l . 

B( I) zH (I )*R*B( I+l) 

21 C( I) z8 (H + R*C( 1+U 

IF ( A BS.( e, (.1.) / HJ_l_)_)r..EL.50,.t_5„(lt_2i 

24 IF (C (2 J) 23 .22. 23 

22 RzR+1, 

__ __G0T0 30 .. 

2 3 RzR-B( 1) /C (2 ) 

30 D037I1ZI.NP 

IzNC-Il 

B( I) ZH (I )-P* 6( I+'l) -a -»B( 1+2 » 



31 

32 
.33 


54 


.75 


77 



37 C( I ) rB ( U-P*C( 1+ 1) -0*C ( 1 + 2 J 
IF (H (2 ) J 32 »3 1. 32 
IF (ASStB(2)/H( U )-E:) 33t33f 34 
IF(A3S(3(2)/H(2))-E)33»3 3*34 
IF ( ASS ( 3 ( 1 ) /H( 1) )-E) 70 »7 0i 34 

34 CBAR rC (2 )-8( 2) 

D-C ( 3) *+ 2-C8 AR *C (4 ) 

IF (D J36. 35 »36 

35 P::P-2.. 

Qz 0* (0 +1 . ) 

■GOTOqg 

36 PzP+ (3 12 ) + C( 3) -8 a J *C« 4) J/D 
0z8+(-3(2)*CBAR + B<i)*rf3n /n 

49 CONTINUE 
E= E * 10 . 

G0T020 

„5D NCzNC-1 — 

VtNCIzO, 

IF ( IRE V) 51 fS2f 52 
.SX.U(NC)rl../f?_., 

G0T053 
52 U(NC)=R 
_S3_D054lz: 1 ».NC-- 
H( I) z8 tI+1 ) 

CONTINUE 

G0T04 - 

70 NCrNC-2 

IF (IREV) 71 t72»72 

-JJ-- GP.rL./.Q 

PPrP/l 0*2.0) 

G0T073 


7 2 GPrO 

PPzP/2 .□ 

7 3 F= (PP) **2-0P 
IF (F )74,.75 t75_ 

74 U(NC + l>z-PP 
U( NOz-PP 
V(NC+1 JzSORT (-( 

V( NC )z-V (NC+11 
G0T076 

U( NC + 1 )z -( PP/ABS (PP) )* ( ABS (P P) +SQRT ( n ) 
V( NC + 1 J = 0. 

U( NC ) = 0P /U (NC+ 1) 

V( NC IrO. 

’‘7 6 D0 77IZ1.NC ' 

H( I) =3 tI+2 ) 

CONTINUE __ 

GO TO 4 

100 RETURN 
FNO 
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Appendix 5 - STVFDBK 
Discussion 

A discussion of this program is given in reference 6. The input 
is formed from the state equations: x = Ax + bu, y « cx where x is the 

state vector,, u is the input and y is the output. The major use of this 
program is to convert the state variable form of the plant to the equi- 
valent open loop transfer functions. 


Input Format 


Card No, 

Column 

Description 

Format 

j 

1 

1-20 

Identification of problem 

4A5 


.21-22 

N = right order of A matrix, an 
integer right justified in the field 

12 

1 

1 

2 

1-10 

~ first element of A matrix 

8E10.0 


11-20 

^12 

• 

• 

• 

8E10.0 

• 

• 

• 

•3 - 

o 

CM 

8E10.0 


11-20 

r 

^22 

• 

• 

• 


n+2 

1 

O 

•^1 

8E10.0 


11-20 

^2 

8E10.0 - 


- 

• 

• 
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Input Format (cont'd) 


Card No. 

Column 

Description 

Format 

n+3 

I-IO 

?1 

8E10.0 


11-20 


8E10.0 

« 


etc. 


• 

• 

n+4 

1-80 

Blank card 



Notes ; 

As many problem may be run as desired by repeating the above set 
of cards for each problem. 
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PROGRAM STVROGK (INPUTt OUTPUT, TAPF<4 z INPUT) 

C OPEN LOOP ONLY 

DIMENSION A ( 1 3» 10 I » Rf 10 ) f H( 10 ) t C ( 1 0 ) » N A ME ( q ) , A A ( 1 0 f 10 ) t C ( 10 ) 

01 ME NS ION 0(lC)»P;i0Tl0)tCC(10)tHH(10).PIM(lG.lD)*U(10)»V{10) 

2000 FOPMAT(qSHO 4 +**«♦♦«•*♦♦♦»***♦♦♦«♦**♦**»*♦*♦♦••+•**»* ) 

2001 FORMAT (2A10»I2) 

2002 format (8E10, 0) 

2003 Format i ihi, sx ,i ^hprobleh ide-nt. ,5x»2jmo) 

2004 format { IHO »5X, 12HTHE A MATRIX/) ' 

20C5 format (6E2D. 7) 

2006 format (1H0»5X, 12HTHE 8 MATRIX/) 

2007 format ( IHO ,5X, 41HTHE CLOSED-LOOP CH A RA CT ER IS T I C PO LY NO HI AL / ) 

2008 FORM AT { IHO ,5 X, 25HTHE FEEDBACK COEFFICIENTS/) 

2009 format ( I’HO i5X» lOHTHE GAIN r CIS, 7) 

2010 FORM AT ( IHO ,5 X, 12HTHE C M AT RI X, 5X , 5H* /j^ 

2D11 FORMAT (IHO fSX, 24HOENOMINATOR COEFFICIENTS/)’' 

2012 FORMAT dHOtSX, 22HNUMERAT0R COEFFICIENTS/) 

2013 FORMAT (IHO r5Xt 23HTHE NUMERATOR OF H- EO UI VA L^NT /) 

2014 FORMAT (IHQ ,5X» 22H0PEN-L00P CALCULATIONS)”"' 

2015 FORM AT ( IHO ,5Xi 2SHMAX IMUH NORMALIZED ERROR = E10.2/) 

_2016_FORMAT (I 1) _ 

2 017 FO RM AT( i HQ , 5 X ,*24 'h'c L 6 SE D~- LO OP CaCcu'l ATI 0 NS ) 

2018 FORMAT { IHD tSXi 6HKEY - 1 1 »3 X, 5H ** ** * ) 

2019 FORM AT ( IHO *5Xt Z3HPLANT IS U NCON TROLL_A BLE5 X » 1C H* ******* *,♦) 

2020 FORK AT (IHO r4Xt’ 35HPLANT "l S“ NUKERKIALL'y UNCO'nTR'o LL AB Le'kd'x »' 

1 16HMAX. DEVIATION z E 1 0 .2 t5 X, 10 H« ♦**♦♦**** ) 

_„2 02 1 . ro RM AT ( IHO , 5 X,14HTHF RO OTS ARE ,13X>9HREAL P AR T, 1 D X t 14 HI M AG IN A ftY„ 
■ lART) 

2022 F0RMAT(25X, 2E20,7) 

2 0.2.3 FORMAT! IHO.) 

C READ INPUT DATA 

1 READ 2001, ( NAME (I ) , 1= 1, 2) ,N 

IF C E Of, » 4.) 1 .0 f„80_0_ 

800 PRINT 2003 »{NAHE(I). Izl,2‘ 

PRINT 2004 

.DO 2 Izl,N : 

READ 2002, ( A (I tJ 1, Jz i, N) 

2 PRINT 2C05 » ( A( I » J) f Jzl »N ) 

PRINT. 20 06 

READ 2 00 2, (B (I ) , Ir.l,N) 

PRINT 2005,(5(1) ,1=1 ,N ) 

_C=_ CHECK. CO NLRO LL.AB IL.1T_Y_ 

DO 7 I=1,N 
AA ( I ,1 ) = B( I) 

J7_ .QONTJNUE. 

DO 8 I=2,N 
L=I-1 

DO a J=1,N 

AA (J ,I ) = 0. 

DO 8 K=1,N 

AA {.J ,I )=.AA ( J ,IJl:Jl(.J».K)JA.A(.K,.U- 

8 CONTINUE 

CONTR=DET{ AA ,N > 

IF(CONTR) 3,4,3. 

4 PRINT 2019 

GO TO 9 

„C NO.TE_USE OF _DUMK.Y. J\R.6U.MENT_C.AL.LELn_ 

3 CALL SIHEQ (AA,C,N,PIN, C) 

DO 43 1=1, N 
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00 t}3 JrltN 
P( It J) -D . 

00 4 3 Kr It N 

P(ItJ)rPtItJ)+AA(ItK?*PIN(Kt.JJ. 
43 CONTINUE 

ERRORzO* 

DO 4 4 in, N 
DO 44 J-ltN 
IF(I-J) 45t46t45 
46, ERR = A3S( P( I. J> -1.0 I . 

GO TO 44 

45 ERR r ABS( P( It J) ) 

4 4 . ERROR - AMAXK ERR t ERROR,). 

IF (ERROR -1 .E-S ) 9t47t47 
47 PRINT 2020tERR0R 

C „ OPEN -LOOP CAL.CUL A.TI.ON.S 

9 PRINT 20 00 
PRINT 2014 
NN -N+1 
PRINT 2011 
CALL CHREOT A tN tD ) . 

.PR I N T 20 05 T ( D.t I ) _t 1 = l_i N N )_ 

CALL PROOT (N rD tU tV T+ 1) 

PRINT 2021 

PRINT 2022 r(U{,I)..V.(I )..f„I=2,_NL 

0011 m rN 
P( 1 1 N) rB (I ) 

-11 CO N I IN UE 

DO 12JJ-2 tN 

0012 m tN 

Jr N- JJ + 1 

Kr J>1 

P{ It J) rO (K )*B( I) 

D 0 1 2- L.r I. , N 

P( It J) rP (I ,J) + A( I.L) *P (L .K > 

12 CONTINUE 

-7-2 READ 2DQ2f. _LCt 

DO 70 I rl t N 
IF(C(D) 7i#70»71 

-70 CONTINUE 

GO TO 1 

71 PRINT 2023 

PRINT 2010 

PRINT 2005 t (C(I)tlrltN) . 

49 0013 1-1 tN 
_-CC (IlrO. 

0013 Jrl ,N 

CC ( I )rP( J, I) *C (J ) + CC (I ) 

-13 CONTINUE 

DO 100 IrltN 
HrNN-I 

JF.(CC(M) )l-01-tl 00 tl.Q) 

100 CONTINUE 

101 PRINt 2012 

..PR IN T -20 05 t ( CC.C II t I.r.1 1 M ) 

MrM- 1 

IF(M) 105fl05tl03 

103-CALL- PRO.OT-(HtCC» Ut-Vt.+ia.- 

PRINT 2021 

PRINT 2022 rfUtI) tV(I )f Izlt M) 



105 


GO JO 12 

10 STOP 
END 

SUBROUTINEPROOTl Nt A, U* V* IR ) 

or ME NS ION A(20)»Ut20)tV{Z0)fHt21»fB(2l)»C(21) 
IREV -IR 
NC rN+1 
DOII::! »NC 
H( I) rA (I ) 

1 CONTINUE 

' . P-0. 

OrO. 

Rz D. 

3 IF (H (1 n tj, 2f 
2 NC:NC-1 

V( NC J = Q. 

“U(NC)rD. 

001002I-1»NC 
H( I) zH(I+l ) 
lOOT' CONTINUE 
G0T03 

4 IF tNC-lJ 5f 100f_5__ 

■5“IF (NC-2‘) If S f Y 

G R=-HU)/H12> 

G0T05D 

' l IF (NC-3‘) 9. S t 9” 

8 PrH( 2) /H [3 ) 

9=H(.l) /H (3J 

GOTO 70 ’ 

9 IF ( ABS IH (NC-I) /H (NO >- A8 5( Ht 2) /H (I n >10# 19*19 

-_1 p_. IR E V z- IR_Ey„ 

MzNC/2 
DOll l3 1, M 

NL -NC 

FrH(NL)'~ 

H(NL)rH(I) 

11 H ( I ) rF 

■ . IF (0)13^12 *13” 

12 Pro. 

.60 TO 15.. 

13 PrP7 0 
Orl . /O 

,_15_.IF'(R )lG.».1.9.#i.6__ 

16 Rrl./R 

19 Ez 5.E-10 
B( NC.)rH( NC )„ 

C( NC )zH ( NO ■ 

B(NC + 1 )rO. 

ONC + l ).rO. 

NPzNC-1 

20 0049Jr 1# 1000 

0021 II rl tNP 

IrNC-Il 

B( I) rH ( I ) + R*B( I* 1) 

,21 C( I) rB (n + R*C{I + .l) 

IF ( A8S ( 3 (1 ) / H( 1) >-E) 50 #5 6* 24 
24 IF (C (2 n 23 *22* 23 

_2 2...R=R+ 1* 

G0T030 

23 R=R-B( 1) /C (2 > 



3D 003711 ::1 .NP -[£2, ' • 

■I-NC-Il 

B(I)rH(I)-P*3CI+D-Q>»3«I+2) 

37 C( I) z3 < I )- P*C( 1+ 1 ) -Q*C ( I-»2 ) 

IF ( H (2 ) ) 32 »3 1. 32 

31 IF (ABS (s3 (2 )/H( 1) )-E) 53 t33» 34 ‘ 

32 IF ( ABS ( B C2 ) /H( 2) )-E) 33 »3 3t 34 

33 IF{ ABSCBCl )/H{ 1) )-E> 70t70t 34 

34 C3AR rC ( 2 3- B( 2) 

Or C ( 53 ** 2- CB AR*C (4 3 
IF (0 33G. 35 »36 

35 PrP-2. 

OrO* (Q+1 .) 

G0T04S 

36 (B(2 3*C( 33 -B U }♦€( 4) 3/0 
Gr 0+ (-3(2) *CBAR-»-B( I) *C .(.3,^) /D 

"4 9 CONTIt^lUE 
EzE* 10. 

2 

”*5 0 NCrNC-1 
V( NC 3 = 0. 

IF { I RE V ) 51 f 3 Qi. 

si U(NC)rir/R' 

G0T053 

_52_U(NC3=R 

5 3 D0 54l-"l» Nc'” 

H(I3rB(I+ll 

_54_1_ .CONTINUE 

.G0T04 

70 NC=NC-2 

IF (I REV) 71 »72Li 72 

71 OPrl./Q 
PPrP/( 0*2. G) 

GOTQJ.3. 

72 OPrO 
PPzP /2-0 

.7.3 Fz(PP3**2-OP 

IF (F )74» 75 r75 
74 U(NC+1) = -PP 

U(NC)r-pP 

V( NC + 1 )rSORT (-FJ 
V( NC )r-V (NC+ 1) . 

G0T0 7G .. 1„. 

75 U( NC + 1 )z-{ PP/ABS (PP) 3* (ABS (PP) +SORT(F) 3 

V( NC + 1 )zO. 

U( NC 3 = GP/U (,NC+.1}_ 

V(NC 3zO, 

76 D077lzi»NC 

:„_H ( I ) =8 ( IjtZJ. 

77 CONTINUE 

60T04 

.. iOLD._ RETURN 

END 

FUNCTION OET ( A »N 3 

_.C. FUNCTION DET DETERMINES T HE.,.DET.ERHINANT.„OF. ..THE. -MA.TR IX -*A* 

DIMENSION A( lOrlO) »8(lDi 10) 

C SET EQUAL TO *A» BECAUSE IrfE DESTROY »A ♦ IN THE PROCESS 

DO..,l_IKz lt_N 

DO 1 JKrl.N 
B(IKtJK) = A(IK»JK) 



n n 


1 continue: 

NN z N-i 
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D z 1. 0 

C IF Nzl THEN BYPASS PIVOT PROCEDURE. .00 DIRECTLY TO CALCULATION 

IF (NN) 69 ,89f 3 ’ - • - ^ 

C START PIVOT SEARCH PROCEDURE 

3 00 100 Lrl .NN 
LL z L-* 1 
AHAX z A (L »L ) 

IH z L 
•JM z L 
DO 15 IzL.N 
DO 15 JzL.N 

IF ( ANAX -ABS (A II »J n >’ 10 *1 5t 15 
10 IM z I 
JM z J 

.AM AX z ABS (A «I .J ) ) 

15 CONTINUE 

found pivot at row ’ IM ♦ ANO_AT COLUMN *JH* _ __ 

IF all REHAININO TERMS "lN'‘THE""MATR ix"A re* zero. sYt ~DET zO AND 
IF tAMAX) 70.70. RETURN, 

C NOW CHANGE ROWS AND rniiiMM^: tp Mrrr<:^ADv 

“14 ' ‘‘“'IF (IM-L)' IsTzO.'i'g 

16 DO 17 Jzl.N 
T z A( IM.J) 

A( IM tJ 1 z A(L, jT 

17 A(L.J) -z T 

. _D z -0. 

20 ~ 1 f'( J M-’l ) ' 2 1 » 25 . 2 . 

21 DO'22 Irl.N 
T z A{.I..JM> 

At I» JM ) z A( 1. L)- 

22 At I.L) z T 

p z. -p 

C PIVOT NOV AT A (L rL r. DIVIDE ROW *L< BY AtL*L) 

25 DO 30 KlzLL.N 

AtL.Kl) r AtL.Kl )/A(L.L) 

30 CONTINUE 

C NOW PRODUCE ZEROES BELOW MAIN DI AG ON AL »A t J tL lz 0, 

DO SO JzLL»,N._ 

DO 50 KzLL.N 

AtJ.K) z A ( J fK )- At J. U *A (L »K 1 

_50„ ^CONTINUE.™ 

100 CONTINUE 

C MULTIPLY MAIN DIAGONAL ELEMENTS TO GET VALUE OF THE DETERMINANT 

63_PQ. 2 00.. I.zl_.N„ 

D z D*At i.I) 

200 CONTINUE 

DET. z 0 . 

C NOV RESTORE. THE VALUES OF THE »A* MATRIX 

DO 2 iKzl.N 

„P0 2 JKzlr.N... 

AtlK.JK) z BtlK.JK) 

2 CONTINUE 

.RETURN 
70 DETzO, 

DO ^ JKzl.N 

DO q lKz-l.,M.. 

AtlK.JK) z B(IK.JK) 

4 CONTINUE 



l62 , " 

RETURN • ■■ 

' END 

SUBROUTINE SIMEO ( A t XO OT »K C » A I NV t X ) 

DIMENSION Al 10»1D) »B«10flO)tXDOT{10) »X (10) »AINV(10»10), . 

DOl 1-1. KC 
DOl Jrl.KC 
AINV ( I . J ) = □. 

1 8 Cl . J) rA ( I . J ) 

002 I-I.KC 
AINV (I .1 )r 1, 

2 X( I) ::X0dT( 1 ) 

D03 1:1. KC 
COMP =0 » 

Kzl 

IF ( ABSla (K »I )) -ABS (COMP) )5 »5 .14 

H COMPrB (K .1 ) 

NzK 

5 K=K+1 

• IF(K-KC) 6.6.7 

't 'IF (0 (M .1 1 ) 8. siVa 

8 IF(N-I)51.i2.3 
_9_DOiO Mz:r.KC 
" 'TEMp'zBd .M) 

8 ( I. M5 =9 (N »M ) 

B( N, M) zTEMP 
TEMPrAlNV( I.H) 

AINV C I .M )z AINV (N »M ) 

.1 D_.A I N V. ( N . M ) :_IEH.P 

TEHPzX (I ) 

X ( I ) rX ( N ) 

X(N)=TE-MP ... 

12 X( I) zX ( I )/8( 1. I) 

TEMP r B (I .1 ) 

.DO 13_ M-::l,.iCC^ 

AlNVd.M) r AINV (I »M )/TEMP 

13 B(I.M) r Bd.M)/TEMP 

DO 16 Jrl.KC 

IF ( J-I )1A. 16 .14' 

14 IF (8 (J .1 ) ) 15 .16. 15 
-1.5~-XC J) :iX ( J B(.J».IJ *X dJ 

-TEMP re (J .1 ) 

0017 Nrl »KC 

. AINV ( J »N I-AINV-IJ-.N )- TE.MP-+A JN V-t-I*-N) 

17 B( J. N) r8 ( J .N )- TEMP*B (I .N ) 

16 CONTINUE 
_..3_(:0N.TINU£__ 

RETURN 

51 PRINT52. I.KC 

32-F0RHAT( .18HQ.LRR0R .I.N„ CO LU MN-.I 2 .. ZX 3 HOX_M A TR I X . » .5X ..3 HK.C^.I2 // ) 
RETURN 
END 

SU B R OU II NE.- CHRE Q-( A., N.. C-) 

DIMENSION Jdl ) . C( 11 ).3 t 10 .1 0) .A ( 1 Q. 10 ) » 01 30D) 

NN rN+1 

DO 20. I rl »NN^ 

C ( I ) rO . 

CONT INUE 

C(.NN 1. -r-_J.. 

0014 M rl »N 
KrO 



Lzl 

JU) rl 
GO TO 2 

1 J ( L) r J ( L) + 1 

2 3f5»50 

3 MMrM-1 

DO ^ T=L»MM 
IIzI+l 

4 J( II >zJt I> +1 

5 CALt: FORM ( J ,N1 ,A »B ) 

K=K+ 1 

0(K) zOET (3 »M J 
DO 6 Izl,M 
LrM- I-t- 1 

IF ( J (LI- (N-M4L n 1*6.50 

6 CONTINUE 

Hi z N-M+l' 

DOm Izl.K 

C(M1 )zC{ Ml )*0<1 »*<-!.) •♦M 

CONTINUE 

RETURN 

50 PRINT 2000 
RETURN 

2000 FORMAT ( IHO. 5X . 1 4HERR0R IN CHREOJ 
END 

SUBROUTINE FORM IJ, H, aVsV 
DIMENSION A ( 10.10) .8 (1 Di 10 ) t J( 11 J 
DO 1 Iz 1, H _ ' 

DOl' KzI.'m 
NR zJ (I ) 

_NCzJ{K)_. 

1 I.K) zA (NR.NC r 
RETURN 
-END 
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Appendix 6 - B0DE4 


Discussion 


This program is used for providing phase and log magnitude versus 
frequency plots for the resulting open and closed loop transfer functions. 
The author is Indebted to Dr, L. P. Huelsman for this program. 


Input Format 


Card No. 

Co lumn 

Description 

Format 


1-80 

hTR =* Title of problem 

80A1 

2 

3 

NC = No. of functions -to be plotted 

13 

3 

1-3 

NP = No. of values of frequency in 
the plot (1-100). 

13 


4-6 

NPD = No. of frequency points/decade 

13 


7-9 

KHZ = of starting frequency in 

Hz 

13 


10-12 

! NSCAL » Max. ordinate of magnitude 
on frequency plot in db. times SCAIM 
j (minimum value is 100 SCALM) db 

lower ) . ! 

13 


13-15 1 

KPHS *= Indicator for phase plot, 0 
if no plot desired, 1 if plot is 
desired 

13 


16-18 

MSGAL =■ Max. ordinate for- phase 
plot in degrees times SCALP 
(minimum ordinate (100 SCALP) 
lower) 

13 


21-30 

SCALM = Scale factor for magnitude 
data (1.0 if no scaling desired) 

ElO.O 


31-40 

SCALP “ Scale factor for phase data 
(1.0 if no scaling desired) 

ElO.O 
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Input Format (cont'd) 


Card No. 

Column 

Description 

Format 

4 

1-3 

ND = Max. degree of numerator and 
denominator polynomials (1-19) 

1 

13 

5 

1-80 

A = Numerator coefficients for 
first function (in ascending order, 
use 10 _ columns for each coeffi- 
cient) 

8E10.0 

6 

1-80 

B = Denominator coefficients for 
first function (in ascending order, 
use 10 columns for each coeffi- 
cient) 

8E10.0 


Notes : 

1 . 

2 . 


3. 

4.. 

5. 

6 . 


Examples 

If the function varies from +10 to -10 db then choose SGALM = 5 ^ 


The maximum magnitude in db which will be plotted is NSCAL/SCALM; 
however, the maximum value indicated on the plot will be NSCAL. 
The minimum magnitude in db which will be plotted is (NSCAL- 100/ 
SCALM; however, the minimum value indicated on the plot is 
NSCAL- 100. 

The above comments also apply to the phase. 

The actual values of the variables are listed in the printout 
for every frequency, 

-rv, ^ . T, . C , nNHZ ^ - _NHZ + NP/ND 

The frequency range in Hz is from 10 to 10 

For each function .to be plotted cards 4 through 6 must be 

repeated. As many cards as necessary may be used to specify 

A and B when ND 8 . 


and NSCAL =50 since 
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NSCAL/SCALM = 10 
(NSCAL-100)/SCALM = -10 

If the minimum frequency is . 01 Hz and the maximum frequency of 
the plot is 100 Hz with 20 points per decade then the total number of 
points is 

NP = decades (points/de cade) = (log 10 (max. freq.) - NHZ) NPD 
= 4 • 20 = 80 



er\ 


PROGRAM B0DE4 ( I NP UT 1 0 UT PU T f TAPE 4-lNPUT) 
DIMENSION LTR(8D) 

DIMENSION Y{5»1D0) 


DIMENSION Xt 10) 

1 PHSOrO. 

READ 100 » LTR 
IF {EOF, 4) 903t2 
2"' •' PRINT 101, LTR - 
READ 105, NC 

PRINT 110, NC _ — 

"lOZ READ 105,NP,NPD,NHZf NSCAL ,NPHS,MSCAL ,SCALM, SCALP 
PRINT 112f NP ,NPD ,NHZ 
DO 140 K=1,NC 
"'PRINT 115, K 

READ 77, (X(I)f 1=1,5) 

77 FORMAT (8F10. 4) 

■ ■'■"PRINT 78»{X{I) .1=1,5)” 

78 FORMAT{* X{D* / lX,iOF12.3) 

NXrl 

"PRINT 201 ” 

00 140 J =1 ,NP 
FR0zXL4(NX,NHZ,NPD) 

RADzG. 28 Z1 e^^FRQ 

CALL MAG77 (RAD.FM, X, PH5D ) 

137 FML=20.^ AL0610(FM) 

' y{ K, J) iFML’*S‘CAL'M 

140 PRINT 20 2, J, FRO , RAO, FH ,FML »P HS D 

PRINT 145, SCALM.LTR 
CALL PL0T4 (Y',NC,NP,NSC'aL) 

GO TO 1 
999 STOP 

iob'>ORMAT' (80A1>“ 

101 FORMAT t IHl, BOA 1) 

105 FORMAT ( Gl 3, 2X »2E1 0. D) . _ / 

110 FORM AT”^'iHO, 16HNUMBER OF CASES=,j.t# 

112 FORMAT {1H0,I3»7H POINTS, 11 0, 7H /DEC AD E, 5X , 1 8H ST AR T I NG FROM 10**, 
113). 

Tl'S F'ORM'aT {1HD,/12H C A S E NU MB ER , 12 ) 

145 FORMAT ( IH 1 , 14HM A6NITUDE DATA ,14H SCALE FACTOR^, E 10 .3 , 1 X * 80 A1 / ! 

201 FO R M AT ( IH 0 , 25 X , 5H HE RT Z 13_X ,_7^AJ3/SEC« 17 X, 3HMAG , 

1/) 

202 FORMAT ■{ IX ,15, 5X ,5E20. B) 

END : 

SUBROUTINE P LO T4 { Y , M ,N F , NS ) 

C -SUBROUTINE FOR PLOTTING 5 X 100 INPUT ARRAY (FORTRAN 4> 

. DIMENSION .Y (.5,10 0) , .LINE (101 ), L( 111., JL (5.) 

DATA ( JL (I ) , 1= 1, 5) / I HA ,1 HB ,1 HC ,1 HO »1HE /, JN » JP, JI f JBL ANK, JZ / 

1 1H-, 1H + , IHI, IH , IHS/ 

_D0 9 9 .1= If 101 

LINE (I )= JBLANX 
. 99 CONTINUE 

N=D 
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C PRINT ordinate SCALE 
DO 101 I-l.ll 
L ( I) C!* I- 11 D+NS 
101 CONTINUE 

PRINT lD5t tL (I ) .Izlf 11 J 
105 format ( 3X » 1 1 ( m »RX) tCHY (I tl ) > 
GO TO 115 

110 IF {N/10-(N-D /ID) 125»125»115 
Q._ ...CONSTRUCT ORDINATE GRAPH LINE 
115 NOrO 

DO 120 IzlflO 


ND-ND+ 1 
LINE (ND) 

DO 120 J=1 »S 
NOrNO+ 1 

120 LINE (NO) rJN 
LINE (101 )rJP 

. IF (N) 135.tl21tl35. 

121 PRINT 170»N»LTNE 
GO TO 185 

JC ^CONSTRUCT 1 „L I N.E 0F_. ABSCI SSA G RAPH LTNFS 

125 DO 1 30 1=1 »101 »10 
LINE (I )=JI 

no. CONTINUE 

C CHANGE NUMERICAL DATa"t 6 Le'tTeRS 
135 DO 160 Irl tM 

XNSzNS _ 

JAzY (I »n')+ fbl. <t3 99S- XNS • 

IF (JA-101) I40f 155* 145 

140 . IF (JA) 150f IS.Of 155„ 

145 LINEdODzJZ 
GO TO 16D 

.15.Q„LINE (1 ) = JZ 

GO TO 160 
155 LINE (JA) zJL(I) 

160 CONTINUE . _ 

C’"' ‘ PRINT LINE of’ data""' 

IF ( N/ 10-{ N- 1) /lO) 175»175»165 

165.. PRINT 17D»N»LINE.Y.(1.»NJ 

170 FORMAT ( IX . 1 4» 10 lA 1 1 IX » E 12 .5 ) 

GO TO 185 

-_17 5.PRINT laOfLINE.Y.Q fN )• . . 

180 format { 5X pIDIAI tlXf E12. 5) 

C SET LINE VARIABLES TO ZERO 

18 5. DO, 130..J.^nJl'01_ 

LINE (I >= J8LANK 
190 CONTINUE 
. .. 195 .NzN+1 

IF (N-NF) 110*11D»200 
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20D RETURN 
END 

FUNCTION XL4 {N»LOiLP> 

AN rN 

ALO-LO 

ALPzLP 

AA zAN/ ALP+ ALO 
IF lAA) lOOtlDDtlOS - 
lOD XL4zl, /ID.** (- AA ) 

50 TO iin 
106 XL4zlO.**AA . 

110 NzN+1 
RETURN 

.. END - - 

SUBROUTINE M AG77 IS »T AMP t Xf PHASE) 

DIMENSION XtlO) 

. COMPLEX GN »GDt ON »DOf PN tCN. CO .T 

GN zCMPLX (79. e» S* 17 2. 0) *CHP LX (7$. S» S- 17 2. 0) * 38Q.0 ♦ 

1 CMPLX (176.D.S+331.D )• * CMPL X ( 17 6. D* S- 33 i. 0) * 

.2 CMPLX (80. 6f S*57S. 0) * CM PL X ( BO .8 » S-S75-.0 J ♦ 

3 CMPLX (82. 8. S*345,0) * CMP LX (3 2. 6» S- 94 5. 0) 

GD zCMPLX ( 46. 0* S) *CMPLX (4 4. 5» S+282. D) *CMPLX (4 4. 5» S-28 2. 0) * 


1 CMPLX (84. 4* S+478. 0) * C MP LX ( 8 4 . 4, S- 47 8. 0) * 

2 CNPLX ( 1 30.5 .S+740.0 ) * CM PL X ( 130, 5* S- 74 0. 0) * 

3 CMPLX (191 .5 .5+1082. 0) * C MP LX ( 1 91 . 5 »S -1 09 2 . 0) 

0Nz5 .1 5E + 19*CMPLX( 10 30.0 »S ) . ... 

DOzCMPLX (650.D»S >*CMPLX( 1260 .0 .S+588 ,D)*CMPLX( 12 60 .0 #S -5 88 .0 ) 

1 *CMPLX(2fl.O.S+147 0.0)*CMPLX(261.0»S -1470.0) 

2 ^CMPLX (391 G. 0* 5+ 18 10 .0 )*CMPLX ( 39 10 .0 »S-1810. 0 

PNr CMPLX( 10 10.0 fS ) / 1010,0 

CNzl .3 35*CMPLX (1 42.. S+ 24 8. )*CMPLX ( 142. .5-248 ,) 

CDzCMPLX (0 .. S) * CMP LX (5 00.. S) 

Tz (GN*DD*CD> /( 6D*D0* CD + GN* 0N*CN*PN ) 

TAMPzCABS( T) 

AzRE AL (T ) 

BrAlMAG( T) 

PHASEzATAN2( B. A) *5 7. 3 
RETURN 

CNzX (3 )*CMPLX( X(^l) ,S) ♦CMP LX (X (4 ) . S) 

CD rCMPLX (X (2 ) . S) *CMPLX ( X(5)»5) *CMPLX(0,.S) 

Tz (GN/60)* (ON/DD)* (CN/CD)* (PN) 

END 
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Appendix 7 - Case Study I 
Second Order Plant 


Discussion 

A second order plant has been selected so that some insight into 
the effect of compensators can be obtained before the computer programs 
are used. Suppose that the plant transfer function is 


6(s) 


100 

s(s+b) 


(6-1) 


and the desired transfer function is a Butterworth type response for 
maximum bandwidth^ 


T(s) 


100 


s +^* lOs+100 


( 6 - 2 ) 


The feedback compensator will be of the form 

K^+H^s 


H(s) - 


Pj^+S 


(6-3) 


where arbitrarily since the niunerator and denominator may be multi- 

plied by .any constant. For zero steady state error at low frequencies 


H^ = Pl 


The transfer function of the series compensator is of the form 


G(s) « 


P^+s 


CDj^+CD^s 


(6-4) 


Prom the relation between the closed loop transfer function and the open 


loop transfer function. 
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X(s) = 


G(sl_C(s) 

l+G(s)C(s)H(s) 


(6-5) 


100 


100(P^+s) 


lOs+100 s(s+b)(CDj^+CD2s)+100(H^+H2S) 


Cross multiplying^ and equating coefficients result in the 
equations, 

3 

s 1 = CD, 


s /r 10+P^ = CDj^+1 CD 2 


lOCH-/? 10 P, 


bCD^+100H2 


Pi = 


Solving for CD2, above equations become 


CD, 


CD, = 


\/T 10+Pj^-b 


H2 = i+CPj^-bxIa'/io - b/ioo) 


H, = P, 


Thus the compensators are specified except for b, and Pj^. 


Sensitivity 

The sensitivity of the closed loop system to variations in plant 

gain is 


S K ^ _T£sI_ 

T G<s)C(s) 2 


s(sd-b) 

10s+100)(Pj^+s) 


(6-6) 
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The low frequency sensitivity is 

b.CD -s 

^ K ^ lOOPj^ 

which is minimum for CD^ = 0 which implies 

CD^ = Jl KH-P^-b = 0 

Pj^ = b-^ 10 


(6-7) 


(6-8) 


Unless P^ is positive the feedback compensator will be unstable and the 
overall system will be infinitely sensitive to variations in Pj^. If Pj^ 
is positive, then (6-8) is a restriction on the plant time constant. The 
only alternative is to make P^^ sufficiently large to keep the sensitivity 
low. cf. eqn. 6~7. It is interesting to compare the sensitivity of this 
system with the Heq compensator and the Guilleman-Truxal compensator. 


Heq Compensator 


Choose Heq to be Heq(s) = l+k 2 S. Then 


T(s) - 


G(s) 


l+G(s)Heq(s) 


100 ^ 100 

s^+10/2 s+100 s(s+b)+ 100 (l+k 2 s) 


Equating coefficients. 


lO-/? b+lOOk^ 

lOv^-b 
2 100 


S 


T 

K 


a 


G(s) 


s (s+b) 

s^+loJ? s+100 
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Gull leman -T r uxa 1 


Choose 


C(s) 


s+b 

s+c 


100 _ 100 

sVloJ? s+100 s(s+c)+100 


Cross multiplying and equating coefficients, 

c = 10/? 


ThuSp the sensitivity is 


3 T ^ s(s-H0/?) 

^ (s^+/? lOs+100) 

The above sensitivities are plotted in Fig. 6-1 for the case of 
where b = 20, and CDj^ = 0. 

If b = 10, then for positive, the sensitivity is minimized 
if is made large. CD^^/Pj^ is nearly one if P^^ « 40. The three sensi- 
tivities are shown below, 

MODOBS 

„ T _ s(s+10)(s-b44.14) 

^ K 2 

(s +10 1^+100) (s+40) 


s 


T 

K 


SS 


s(s+10) 

(s^+10v5“s+100) 




Frequency (Hz.) 


Figure 6.-1 Sensitivity of High-perfoirmance Plant 
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Gullleman-Truxal 

S o 5(6+10^2) 

^ (s\l0^s+100) 


Notice that for b = 20 the MODOBS system is significantly better 
at all frequencies below iHz. The Heq compensation is both unrealizable 
and largest in sensitivity for this case. When the forward gain is small 
then no compensation will help the sensitivity. Thus all the proposed 
systems have a sensitivity near one for frequencies above iHz. 

For b = 10, the low frequency plant gain must be greater to 
achieve the same T(s). The Guilleman-Truxal compensation decreases the 
forward gain proportionately, so that the sensitivity is unchanged. The 
Heq compensation does not change the forward gain; therefore the sensitivity 
is reduced from the case where b = 20. The MODOBS system also improves 
in sensitivity, but this system will only be as low in sensitivity as the 
Heq method when is infinitely large. Nevertheless, the MODOBS 
compensation has the advantage of being realizable. 


Integral Sensitivity 

For the plant with a pole at -20, the integral sensitivity 
(weighted by 1/s) is 


S 


JL 

2rt 


-joo 



^ s(s+20)(s+Pj^-5.86) 

® (s^+loV?" 6+100) (s+P^) 


2 


ds 


This integral can be evaluated directly from Newton, Gould and Kaiser, 
p. 372. 

The minimum low frequency sensitivity was attained for » 5.86 


which from M0D0BS3 has an integral sensitivity of 0.12. The minimum 
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integral sensitivity is 0.085 for = 8. In this simple case this can be 
determined analytically by plotting S as is varied. 

Final Design 

From the integral sensitivity data, a final design choice is 

b = 20 

Pi = 

The following graphs are shown in Figs. 6-2 to 6-4. 

Fig. 6-2 Plant transfer function 

Fig. 6-3 Sensitivity function 

Fig. 6-4 Loop gain function 

Conclusion 

From the above figures the design demonstrates the dependence of 
the sensitivity on the plant parameters. If the time constants in the 
closed loop transfer function are much greater than those in the plant, 
a sacrifice in the system sensitivity at low frequencies must be made. 

In addition the compensator will have large gains associated with the 
feedback which complicate the design of the compensators and increase 
the effect of the measurement noise. 



2 



Frequency (Hz.) 


Figure 6.2 Open and Closed Loop Transfer Functions 




Frequency (Hzo) 


Figure 6*4 Loop Gain Function of Second Order System 


VO 
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Appendix 8 - Case Study II 
Third Order Plant 


Discussion 

A third order plant has been selected to compare the computer 
results with hand calculations. The results agree in all the calculations 
discussed below. The plant, transfer function is chosen to be of the form 


G(s) 


1000. 

s^+GD^s^+GD^s+GD^ 


and the desired closed loop transfer function is again chosen with a 
Butterworth type denominator. 


T(s) 


1000 

s\20s^+200s+1000 


M0D0BS2 

2 

The equations for the coefficients of P(s)“P^+P 2 ®'*'^ 3 ® which give 
minimum low frequency sensitivity were derived in Appendix 2. For N‘=3, 


P3 » 1 

» GD3 - 20 

Pj, <= GD^ - 200 - 20<P2> 

These equations are satisfied if 


GD, = 20 
S 

GD^ =*200 

Low frequency sensitivity will be improved if the plant is ’’Type 1” or 
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equivalently, has a pole at the origin; thus 

GD^ = 0 

The poles of the plant are the roots of 

s(s^+20e+200) = 0 
s = -10 + jlO 

For the plant to have real roots then, 

A 

GD^ - 4GD2 > 0 

and in order for and to be positive, 

GD^ > 20 

GD^ > 200 + 20?2 

These requirements can be met if 

GDg > 68.28 

20(60^ - 10) ^ GD^ <00^/4 
GD^ = 0 

If GD^ is near 68 then must be near zero. Under these circum- 
stances the zero of CD(s) cancels the pole at a low frequency, and the 
low frequency sensitivity suffers. Suppose that 

GD^ =100 




Figure 8.1 Sensitivity o£ Third Order System 
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Then 

1800 < GD 2500 
2 


P3 

= 1 


P2 

= GD 3 

- 20 « 80 

PI 

= gd^ 

- 1800 


To keep the low frequency pole of P(s) as large as possible it is best to 
make GD2 large; thus 

GD3 = 100 
GD2 = 2500 
GDI == 0 
P3 = 1 
P2 => 80 
PI ■= 700 


Sensitivifs 


If the plant can be selected as above, the sensitivity becomes 


K G(s)C(s) 


s(s +1005+2500) s 
(s^+20s^+200s+1000)(s^+80s+700) 


For this choice of the plant the system will have minimum low 
frequency sensitivity and a stable feedback compensator. A comparison 
of the sensitivity of this system with only a series compensator which 
achieves the same closed loop response is shown in Fig. 8-1. Except 
for frequencies between 1.41 Hz and 14.1 Hz the sensitivity is Improved 
by using the MODOBS system. 
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Discussion 

Again, the importance of choosing a plant which is consistent 
with the requirements of the overall system is evident. A slow system 
cannot be made to react quickly and also have good sensitivity charac- 
teristics using second order compensators. 
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Appendix 9 

Case III In core Thermionic Reactor 

The plant and the desired transfer function for an in core 

[7I 

thermionic reactor have been discussed by Weaver et al . The charac- 
teristics of the system are given in the above report. This appendix 
addresses itself to the question of realization of the desired transfer 
function for the fourth order model when 'only the output is available. 
The plant transfer function is given by 

p, . (s+.Q326)(s+.lKs+1.79y 200 

“ (s+.0119+j.0l64)(s+.0119-j.0164)(s+.00663)(s+1.85) 

and the desired transfer function is 


Y(s) 200(s+.Q326)(s+.l)(s+1.79) 

R(s) “ 100(s+.0326)(s+.l)(s+1.79).(s+2) 


When these terms are multiplied out the resulting transfer functions are 

. 1. 16708+48. 1228s+384.52s\200s^ 

G(s) = 7 7 y r — r 

5.03585*10 +1.0542*10 +. 0568639s +1. 8804s +s 

Y(s) ^ 1. 16708+48. 1228S+34. 523^+2003^ 

1. 16708+48. 70635s+408.5815s^+392.26s’\l00s^ 


This problem represents a special case' in that the numerator of G(s) and 
Y/R^°^ is the same. Equating the real and desired closed loop transfer 
functions yields 


T(s) 




1+G(s)C(s)H(s) 


and in terms of the numerator and" denominator polynomials. 
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TN(s) ^ GN(s)CN(s)HD(s) 

TD(s) “ GD(s)CD(s)HD<s)+GN(s)CN(s)HN(s) 


Since CN(s) = HD(s) ^ P(s) and TN(s) = GN(s), 


.... 1.. , = g(s) 

TD(s) GD(s)CD(s)+GN(s)HN(s) 


Equating the like terras in the above equation and expressing the 


in vector matrix form, one gets 


GD(1) 

0 

0 

0 


CD(1) 


GN(1) 

0 

0 

0 

GD(2) 

GD(1) 

0 

0 


CD(2) 


GN(2) 

GN(1) 

0 

0 

GD(3) 

GD(2) 

GD(1) 

0 


CD(3) 


GN(3) 

GN(2) 

GN(1) 

0 

GD(4) 

GD(3) 

GD(2) 

GD(1) 


CD(4) 


GN(4) 

GN(3) 

GN(2) 

GN(1) 

GD(5) 

GD(4) 

GD(3) 

GD(2) 


+ 

0 

GN(4) 

GN(3) 

GN(2) 

0 

GD(5) 

GD(4) 

GD(3) 


0 

0 

GN(4) 

GN(3) 

0 

0 

-GD(5) 

GD(4) 


0 

0 

0 

GN(4) 

0 

0 

0 

GD(5) 


0 

0 

0 

0 


TD(1) 0 

TD(2) TD(1) 

TD(3) TD(2) 

TD(4)- TD(3) 
TD(5) TD(4) 

0 ,'TD(5) 

0 0 
0 0 


0 

0 

TD(1) 

TD(2) 

TD(3) 

TD(4) 

TD(5) 

0 


0 

0 

0 

TD(1) 

TD(2) 

TD(3.) 

TD(4) 

TD(5) 


P(l) 

P(2) 

P(3) 

P(4) 


result 

HN(1) 

HN(2) 

HN(3) 

HN(4) 
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Plugging in the appropriate coefficients from the above expressions for 
TD(s), GN(s), and GD(s) the unknowns GD(s) and HN(s) are equal to 


CD(1) 


CD(2) 


1.191L4E+01 

-5.44224E+00 

9.63714E+00 

-1.72292E+01 


P(l) 



1.41485E+02 

-2.12494E+02 

3.91935E+02 

-7.00789E+02 


« 


9.32689E+02 

-1.65161E+03 

Z.96273E+03 

-5.28468E+03 


P(2) 

• 

= 

3.31178E-14 

3.56005E-14 

-1.10191E-13 

l.OOOOOE+01 


HN(1) 


9.99949E-01 

2.34842E-05 

-4.15859E-05 

7. 4346 7E- 05 


P(3) 

• 


4.963L9E-01 

1.00486E+00 

-8.68215E-03 

1.55222E-02 


• 


5.09230E-01 

7.61501E-01 

5.35279E-01 

8.30837E-01 


P(4) 

« 


4.66345E+00 

8.25805E+00 

-1.47636E+01 

2.82909E+01 



HN(4) 


The resulting expression can also be solved for P(s) in terms of C(s)^ that 
iS; 


P(l) 




CD(1) 

P{2) 


1.10028E-01 1.28262E-03 -5.27575E-04 6.47151E-04 


CD(2) 



-3.32179E-01 1.62916E-01 -2.04714E-02 2.61570E-02 



P(3) 


-2.19815E-01 9.04156E-02 -1.09084E-02 1.92750E-01 


CD(3) 



000 l.OOOOOE-01 



P(4) 


— — 


CD(4) 


If one selects CD(1) equal to zero then there will be an integrator in the 
forward path which results in a system with zero low frequency sensitivity. 


The requirements for open loop stability are positive coefficients and 


CD(1) CD(4)<CD(2) CD(3) 
P(l) P(4) < P(2) P(3) 


These requirements can easily be met when CD(1) is zero if 


CD(3) < CD(2) + CD(3) 
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The choice &r CD(s) is arbitrarily made 

2 

GD(s) = s(s +S+1) 

then 

P(s) = 7.227-10"^ + .1686 s + .27226 + .1 

The resulting system has a stable open loop compensator and zero low 
frequency sensitivity and zero steady, state error. These characteristics 
are an improvement over the response of the state variable feedback 
system while both have the same closed loop response. In addition only 
the output state needs to be measured; however, two third order compen- 
sators are required. 



